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Abstract 

The objectives of this research are: (1) to develop and implement a new methodology for 
large eddy simulation of (LES) of high-speed reacting turbulent flows. (2) To develop algebraic 
turbulence closures for statistical description of chemically reacting turbulent flows. We have 
just completed the third year of Phase III of this research. This is the Final Report of our 
activities on this research sponsored by the NASA LaRC under Grant NAG-1-1122. 


Technical Monitor 


Dr. J. Philip Drummond (Hypersonics Air-Breathing Propulsion Branch, NASA LaRC, Mail Stop 
197, Tel: 757-864-2298) is the Technical Monitor of this Grant. 


I Summary of Achievements 

We just completed Year 3 of the Phase III activities on this NASA LaRC sponsored project (Grant 
NAG-1-1122). The total time allotted for this phase is three years; this phase was followed at the 
conclusion of Phase II activities (also for three years). In total we have completed 9 years of LaRC 
supported research. This is the Final Report and provides a summary of our accomplishments in 
Phase III of this research (August 1, 1996 - July 31, 1999). Our work in this phase consists of the 
following two constituents: (1) development of LES methodologies via probability density function 
(PDF) methods and numerical solution of the PDF via Monte Carlo schemes, , and (2) Development 
of algebraic turbulence closures for statistical description chemically reacting turbulent flows. Each 
of these two constituents are discussed below in order: 
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1.1 LES of Turbulent Reacting Flow 


A methodology termed the “filtered density function” (FDF) is developed and implemented for LES 
of chemically reacting turbulent flows [1]. In this methodology, the effects of the unresolved scalar 
fluctuations are taken into account by PDF of subgrid scale (SGS) scalar quantities. A transport 
equation is derived for the FDF in which the effect of chemical reactions appears in a closed form. 
The influences of scalar mixing and convection within the subgrid are modeled. The FDF transport 
equation is solved numerically via a Lagrangian Monte Carlo scheme [2] in which the solutions of 
the equivalent stochastic differential equations (SDEs) are obtained. These solutions preserve the 
Ito-Gikhman nature of the SDEs. The consistency of the FDF approach, the convergence of its 
Monte Carlo solution and the performance of the closures employed in the FDF transport equation 
are assessed by comparisons with results obtained by direct numerical simulation (DNS) and by 
conventional LES procedures in which the first two SGS scalar moments are obtained by a finite 
difference method (LES-FD). These comparative assessments are conducted by implementations of 
all three schemes (FDF, DNS and LES-FD) in a temporally developing mixing layer and a spatially 
developing planar jet under both non-reacting and reacting conditions. In non-reacting flows, the 
Monte Carlo solution of the FDF yields results similar to those via LES-FD. The advantage of 
the FDF is demonstrated by its use in reacting flows. In the absence of a closure for the subgrid 
scalar fluctuations, the LES-FD results are significantly different from those based on DNS. The 
FDF results show a much closer agreement with filtered DNS results. The results of this work are 
published in Refs. [3,4] and are provided in Appendix I and Appendix II. 

We have also developed a methodology, termed the “filtered mass density function” (FMDF) for 
LES of variable density chemically reacting turbulent flows. This methodology is based on the 
extension of the FDF and represents the joint probability density function of the subgrid scale 
(SGS) scalar quantities. The FMDF is obtained by solution of its modeled transport equation. 
In this equation, the effect of chemical reactions appears in a closed form and the influences of 
SGS mixing and convection are modeled. The stochastic differential equations (SDEs) which yield 
statistically equivalent results to that of the FMDF transport equation are derived and are solved 
via a Lagrangian Monte Carlo scheme. The consistency, convergence, and accuracy of the FMDF 
and the Monte Carlo solution of its equivalent SDEs are assessed. In non-reacting flows, it is shown 
that the filtered results via the FMDF agree well with those obtained by LES-FD. The advantage of 
the FMDF is demonstrated in LES of reacting shear flows with nonpremixed reactants. The FMDF 
results are appraised by comparisons with data generated by direct numerical simulation (DNS) and 
with experimental measurements. In the absence of a closure for the SGS scalar correlations, the 
results based on the LES-FD are significantly different from those obtained by DNS. The FMDF 
results show a closer agreement with DNS. These results also agree favorably with laboratory 
data of exothermic reacting turbulent shear flows, and portray several of the features observed 
experimentally. This work is described in detail in Appendix III. This appendix is scheduled to be 
published in Journal of Fluid Mechanics [5]. 
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We have recently extended the FDF methodology to also include the velocity field. Some prelim- 
inary results are provided in Appendix IV. In this portion of our work, the joint-velocity filtered 
density function of the velocity (VFDF) is considered. A transport equation is proposed for the 
VFDF in which the unclosed terms are modeled in a manner analogous to conventional second 
order subgrid scales closures [6]. The modeled VFDF transport equation is solved numerically via 
a Lagrangian Monte Carlo scheme in which the solutions of the equivalent stochastic differential 
equations are obtained. The consistency and the convergence of the simulated results are assessed 
by comparison with results obtained by LES-FD in which the equivalent transport equations of the 
subgrid scale moments are solved. The accuracy and reliability of the velocity FDF are also assessed 
via comparison with DNS and experimental data. The DNS data are those for a three-dimensional 
(3Z?) temporal mixing layer, and experimental data are those of a 3 D turbulent jet. 

1.2 Algebraic Modeling 

In this work, explicit algebraic scalar flux models which are valid for three-dimensional turbulent 
flows are derived from a hierarchy of second-order moment closures. The mathematical procedure is 
based on the Cayley- Hamilton theorem and is an extension of the scheme developed by Taulbee [7]. 
Several closures for the pressure-scalar gradient correlations are considered and explicit algebraic 
relations are provided for the velocity-scalar correlations in both nonreacting and reacting flows. 
In the latter, the role of the Damkohler number is explicitly exhibited in isothermal turbulent flows 
with nonpremixed reactants. The relationship between these closures and traditional models based 
on the linear gradient diffusion approximation is theoretically established. The results of model 
predictions are assessed via comparison with available laboratory data in turbulent jet flows. This 
work is published in Ref. [8], which is included here as Appendix V. 

The extension of the methodology above for high speed flow has also been completed. In this part 
of our work, closure for the compressible portion of the pressure-strain covariance is developed. It 
is shown that, within the context of a pressure-strain closure assumption linear in the Reynolds 
stresses, an expression for the pressure-dilatation can be used to construct a representation for the 
pressure-strain. Additional closures for the unclosed terms in the Favre-Reynolds stress equations 
involving the mean acceleration are also constructed. The closures accommodate compressibility 
corrections depending on the magnitude of the turbulent Mach number, the mean density gradient, 
the mean pressure gradient, the mean dilatation, and, of course, the mean velocity gradients. 
The effects of the compressibility corrections on the Favre-Reynolds stresses are consistent with 
current DNS results. Using the compressible pressure-strain and mean acceleration closures in the 
Favre-Reynolds stress equations an algebraic closure for the Favre-Reynolds stresses is constructed. 
Noteworthy is the fact that, in the absence of mean velocity gradients, the mean density gradient 
produces Favre-Reynolds stresses in accelerating mean flows. Computations of the mixing layer 
using the compressible closures developed are described. Favre-Reynolds stress closure and two- 
equation algebraic models are compared to laboratory data for the mixing layer. Experimental data 
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from diverse laboratories for the Favre- Reynolds stresses appears inconsistent and, as a consequence, 
comparison of the Reynolds stress predictions to the data is not conclusive. Reductions of the kinetic 
energy and the spread rate are consistent with the sizable decreases seen in these classes of flows. 
Appendix VI provides a complete description of this portion of our activities. This Appendix is to 
be published in Physics Fluids in September 1999. 
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A methodology termed the “filtered density function” (FDF) is developed and implemented for 
large eddy simulation (LES) of chemically reacting turbulent flows. In this methodology, the effects 
of the unresolved scalar fluctuations are taken into account by considering the probability density 
function (PDF) of subgrid scale (SGS) scalar quantities. A transport equation is derived for the FDF 
in which the effect of chemical reactions appears in a closed form. The influences of scalar mixing 
and convection within the subgrid are modeled. The FDF transport equation is solved numerically 
via a Lagrangian Monte Carlo scheme in which the solutions of the equivalent stochastic differential 
equations (SDEs) are obtained. These solutions preserve the Ito-Gikhman nature of the SDEs. The 
consistency of the FDF approach, the convergence of its Monte Carlo solution and the performance 
of the closures employed in the FDF transport equation are assessed by comparisons with results 
obtained by direct numerical simulation (DNS) and by conventional procedures in which the 
first two SGS scalar moments are obtained by a finite difference method ( LES -FD) These 
comparative assessments are conducted by implementations of all three schemes (FDF, DNS and 
LES-FD) in a temporally developing mixing layer and a spatially developing planar jet under both 
non-reacting and reacting conditions. In non-reacting flows, the Monte Carlo soludon of the FDF 
yields results similar to those via LES -FD. The advantage of the FDF is demonstrated by its use in 
reaedng flows. In the absence of a closure for the SGS scalar fiuctuadons, the LES-FD results are 
significantly different from those based on DNS. The FDF results show a much closer agreement 
with filtered DNS results. © 1998 American Institute of Physics. [SI 070-663 1(98)01 402-0] 


I. INTRODUCTION 

Over the past 30 years since the early work of 

Smagorinsky, 1 significant efforts have been devoted to large 

eddy simulation (LES) of turbulent flows. 2 * * * * * 12 The most 

prominent model has been the Smagorinsky eddy viscosity 

closure which relates the unknown subgrid scale (SGS) Rey- 
nolds stresses to the local large scale rate of flow strain. 13 
This viscosity is aimed to provide the role of mimicking the 
dissiparive behavior of the unresolved small scales. The ex- 
tensions to “dynamic” models 14,13 have shown some im- 
provements. This is particularly the case in transidonal flow 
simularions where the dynamic evaluations of the empirical 
mode! “constant” result in (somewhat) better predicrions of 
the large scale flow features. 

A survey of combustion literature reveals relatively little 
work in LES of chemically reacting turbulent flows. 7,16 It 
appears that Schumann 17 was one of the first to conduct LES 
of a reacting flow. However, the assumption made in this 
work simply to neglect the contribution of the SGS scalar 
fluctuations to the filtered reaction rate needs to be justified 
for general applications. The importance of such fluctuations 
is well recognized in Reynolds averaged procedures in both 
combustion 18 " 20 and chemical engineering 21 " 24 problems. 
Therefore, it is natural to believe that these fluctuations are 


also important in LES. McMurtry et a/., 25 * 26 Sykes et a/., 27 
Liou et aL, 28 Menon et al ™ Boris et a/., 30 Fureby et ai , 31,32 
Cook et aL 3334 Mathey and Chollet, 35 Branley and Jones 36 
and others provide several means of conducting LES of tur- 
bulent reacting flows. 

Modeling of scalar fluctuations in Reynolds averaged 
methods has been the subject of broad investigations since 
the pioneering work of Toor. 37 An approach which has 
proven particularly useful is based on the probability density 
function (PDF) or the joint PDF of scalar quantities. 38 " 41 The 
systematic approach for determining the PDF is by means of 
solving the transport equation governing its evolution. 42 In 
this equation, the effects of chemical reaction appear in a 
closed form; this constitutes the primary advantage of the 
PDF schemes in comparison to other statistical procedures. 
The use of PDF for LES was suggested by Givi 7 and its first 
application is due to Madnia and Givi. 43 In this work, the 
Pearson family of distributions are assumed to characterize 
PDF of SGS scalars in homogeneous flows under chemical 
equilibrium conditions. This procedure was also used by 
Cook and Riley. 44 The extension of assumed PDF methods 
for LES of non-equilibrium reacting shear flows is reported 
by Frankcl et al. 45 While the generated results are encourag- 
ing, they do reveal the need for more systematic schemes in 
which the transport of the PDF of SGS scalar quantities are 
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considered. Pope 16 introduced the concept of “filtered den- 
sity function’ * (FDF) which is essentially the PDF of SGS 
scalar variables. With a formal mathematical definition of the 
FDF, Pope 16 demonstrates that the effects of chemical reac- 
tion appear in a closed form in the FDF transport, thus mak- 
ing it a viable candidate for LES of chemically reacting 
flows. Gao and O'Brien 46 develop a transport equation for 
the FDF and offer suggestions for modeling of the unclosed 
terms in this equation. 

The objective of the present work is to further demon- 
strate the applicability of the FDF and to provide results 
based on its implementation for LES of chemically reacting 
turbulent flows. Only the FDF of scalar quantities is consid- 
ered; probability treatment of the subgrid velocity fluctua- 
tions is postponed for future work. 


II. FORMULATION 

We consider an incompressible (unit density), isother- 
mal, turbulent reacting flow involving N , species. For the 
m nth^matiral description of this flow, the primary transport 
variables are the velocity vector u t (x,t) (« = 1,23). the pres- 
sure p(x,t), and the species’ mass fractions 4> a (x,t) (a 
= 13, The equations which govern the transport of 
these variables in space (x,) and time (r) are 


0. 

dXi 


Suj 

dt 

dt 


dUjtij _ dp dTjj 


dXi 

SUj<t> c 

dXi 


dxj dXi ' 

■*r 


a) 


( 2 ) 


(3) 


where w 0 (x,r)=o> a (4»(x,r)) denotes the chemical reaction 
term for species a, and [ d>i ,4>i , . . ■ ,0^] denotes the 

«r«iar array. Assuming a Newtonian flow with Fick’s law of 
diffusion, the viscous stress tensor r t j and mass flux J“ are 
represented by 


dttj du 


dXi 


dXi 


d<f> a 

j a =-r— 

J ‘ d Xi ’ 


(4) 


where v is the fluid viscosity and T is the diffusion coeffi- 
cient, T = v/Sc, and Sc is the molecular Schmidt number. 
Large eddy simulation involves the use of the spatial 


_47 


filtering operation 

(f(x,t)) L = j f(x',t)9(x',x)dx’. 


(5) 


where ^denotes the filter function, (f(x,t)) L represents the 
filtered value of the transport variable /(x,r), and f -f 
-(f) L denotes the fluctuations of / from the filtered value. 
We ™n«idw spatially and temporally invariant and localized 
filter functions, thus ?(i' i «)»G(i , -x) with the 
properties, 47 G(x) = G( — x), and / _„G(x)dx = 1. More- 
over, we only consider “positive” filter functions as defined 
by Verm an et al. 4 * for which all the moments fZ w x m G(x)dx 
exist for ms* 0. The application of the filtering operation to 
the transport equations yields 


6Xi 


=0. 


( 6 ) 


d(Uj) L d(Uj)i(Uj) L <?(p)t | T i ;) L dTjj 

dt <?x, dxj dXj dx, 


dt dXi 


d(J?)L *M° 

— 7— + <"«>!. • 

dx t dXi 


( 8 ) 


where Tij=(uiUj)L—(Uj)i(Uj)i and 
— ( Ui ) L ( <l>a)L denote the subgrid stress and the subgrid mass 
flux, respectively. 


III. CLOSURE STRATEGY 

In LES of non-reacting flows the closure problem is as- 
sociated with 3 Tij=(uiUj) L —(ui)L(uj)L M*=(ui<t> a )L 
k reacting flows, an additional model is re- 
quired for . Here, modeling of ** die subject of 
the probability formulation as described in the next section. 
For the former two, we make use of currently available clo- 
sures which are well-established in non-reacting flows. The 
subgrid stress is modeled via 

T ij -(S ij f3)T kl =-2v l (S ii ) L , (9) 

where (S^l is the resolved strain rate tensor and v t is the 
subgrid viscosity. We use two closures to represent this vis- 
cosity. The first is the same as that in the conventional Sma- 
gorinsky closure 3 

Vl =C,A 2 GyKSij)L(Sij)L, (10) 

where is the filter and C s is an empirical constant 
The drawbacks of this closure arc well-recognized. 49,50 In an 
attempt to overcome some of these drawbacks, we also make 
use of a second closure in which the subgrid viscosity is 
determined based on the modified subgnd kinetic energy 

v ( =c t A C vi(«r>t{«r>z-«“r)t>L'«“r>t)ti- < n > 

where u* = u i —‘26 i and $5^ is a reference velocity in the x, 
direction. The subscript L' denotes the filter at the secondary 
level which has a characteristic size (denoted by A c .) larger 
than that of grid level filter. This model is essentially a modi- 
fied version of that proposed by Bardina et of., 51 which uti- 
lize equal sizes for the grid and secondary filters. We refer to 
this as the modified kinetic energy viscosity (MKEV) clo- 
sure. 

A similar model is used for the closure of the subgrid 
mass fluxes 32 


A#r=-r, 


dx t ’ 


( 12 ) 


where r,= v,/Sc,. and Sc, is the subgrid Schmidt number 
and is assumed constant. 
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IV. FILTERED DENSITY FUNCTION (FDF) 

The key point in this formulation is to consider the scalar 
fluctuations of the underlying scalars' array 4>(x,r) in a 
probabilistic manner. For that, we define the “filtered den- 
sity function” (FDF), denoted by P Ll as 16 

PJV-.x,/)- J_ + Je[^,4>(x',r)]C(x'-x) < ix'. (13) 

", 

e[*,*(x,/)]=*'i'-<i»(x,/)]=n 

tf* 1 

(14) 

where S denotes the delta function and ^ denotes the com- 
position domain of the scalar array. The term g[<I>- V(x^)] 
is the “fine-grained” density, 39 * 40 and Eq. (13) implies that 
the FDF is the spatially filtered value of the fine-grained 
density. Thus* Pi gives the density in the composition space 
of the fluid around x weighted by the filter G. With the 
condition of a positive filter kernel, 4 ® P L has all the proper- 
tics of the PDF. 40 

For further developments, it is useful to define the “con- 
ditional filtered value” of the variable Q(x,t) by 


<(2(x,r)|^>^ 


/!:G(x'.Oet^.»I»(x'./)]G(x'-x)dx' 

Pt^ix,/) 


(15) 

where {a\$) L denotes the filtered value of a conditioned on 
fi. Equation (13) implies 

(i) For G(x,/) = c, (G(x,r)|^) t =c, 


(iii) Integra] property: (Q(x,t)) L 




(Q(xj)\H') L P L (V;x,t)dV f 


(18) 


where c is a constant, and £(<I>(x,f))*C(x,r) denotes the 
case where the variable Q can be completely described by 
the compositional variable 4>(x,f). From these properties it 
follows that the filtered value of any function of the scalar 
variables (such as the reaction rate) is obtained by integration 
over the composition space 

<G(*. 0h= l + _~Q(*)P L (*;x,t)dV. (19) 

To develop a transport equation for the FDF, the time- 
derivative of Eq. (13) is considered 

f- 3<t> a (x',t) <?e[¥,*(x',/)] 


(*;»,/) f" S^gjJ 

St J —m St 


dlpa 


J) 


xG(x'-x)dx' 

Sl/l a J-m St 

xe[V,<t »(x',f)]G(x'— x)dx\ (20) 


where the summation convention applies to the species suf- 
fix, a. This combined with Eq. (15) yields 


SP L (V;x,t) 

St 


3 


l3<t> a 

,\ , 1 

[(-=- 

*] P L (V\x,t) 


(21) 


Substituting Eq. (3) into Eq. (21) yields 


SP L (V;x,t) 3 


St 


\ 


V) + 

L 


SXi 


-<«.(«I>)|V> L 


Pi(V;x.f) 


( 22 ) 


in which the convective term can be represented in the form 

3 

sWi 


dXi 


(23) 


The unclosed nature of convection is denoted by the condi- 
tional filtered value of the velocity which is further decom- 
posed into resolved and subgrid scale components. It is use- 
ful to adopt the decomposition 

so that Eq. (21) can be expressed as 

SPl S{ui) L P L 5[(«,l 1 i r >t”<«,)JPt 
St 


(24) 


SXi 


SXi 


(16) 

d 

[(£ 

V) P, 

(17) 

Wa 


h J 


*«.(*)J»J 

3<P a 


(25) 


This is an exact transport equation for the FDF and is similar 
to that presented by Gao and O'Brien. 46 The last term on the 
right hand side of this equation is due to chemical reaction 
and is in a closed form. The second term on the left hand 
side represents the filtered convection of the FDF in physical 
space and is also closed (provided (u,) t is known). The un- 
closed terms are associated with the first term on the right 
hand side denoting the effects of unresolved subgrid scale 
convection, and the second term on the right hand side rep- 
resenting the influence of molecular diffusion. 

The subgrid convective flux is modeled via 




(26) 


The advantage of the decomposition (Eq. (24)) and the sub- 
sequent model (Eq. (26)) is that they yield results similar to 
that in conventional LES for the first moment of the FDF. 
The first moments corresponding to Eqs. (24) and (26) are 

*(*a)L 




Sx t 


(28) 
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‘ rfX l (/) = Z) 1 (X(r) l r)d/+£(X(/),i)rfW i (f), (35) 

where X { is the Lagrangian position of a stochastic particle, 
D { and £ are known as the “drift” and “diffusion” coeffi- 
cients, respectively, and W, denotes the Wiener-Levy 
process. 65 A comparison of the Fokker-Plank equation cor- 
responding to Eq. (35) with the FDF transport equation (32) 
determines the appropriate specification of the coefficients to 
be 


, *(r+r,) 

£-V2(r + r f ). D,m(u,) L +- - ( 36 ) 

dx i 

Thus the SDE which represents the spatial transport of the 
FDF is 


dXi(t) = 


<«,>! + 


fCT + TJ 

dx, 


<fr+[2(f + T ,)] m dW j . 


(37) 


The compositional makeup of the particles evolves simulta- 
neously due to the actions of subgrid mixing and reaction 


d<K 

dt 


-n m (*;-<* a > t )+o, a , 


(38) 


where =^ a (X,(/),f) denotes the scalar value of the par- 
ticle with the Lagrangian position vector X { . 

In the numerical implementation, the FDF is represented 
by an ensemble of Monte Carlo particles, each with a set of 
scalars ^ ) (X < " ) (0-0 and Lagrangian position vector X <n> . 
Numerically, a splitting operation is employed in which the 
transport in the physical and the compositional domains are 
treated separately. The simplest means of simulating Eq. (37) 
is via the Euler-Maruyamma approximation 66 

X-(/ i+1 )=X-(/ i )+D-(r t )Af+£"( t4 )(Ar) 1 ' 2 f-(r t ), 

(39) 

where D*(t k )=Di(X M (t k ),t), E H (t k )=E(X { * ) (t k ),t) and 

is a random variable with the standard Gaussian PDF. 
This formulation preserves the Markovian character of the 
diffusion process 67-69 and facilitates affordable computa- 
tions. Higher order numerical schemes for solving Eq. (37) 
are available, 66 but one must be cautious in using them for 
LES since the diffusion term in Eq. (35) depends on the 
stochastic process X(f). The numerical scheme must pre- 
serve the Ito-Gikhman 70,71 nature of the process. The coeffi- 
cients D, and £ require the input of the filtered mean veloc- 
ity and the diffusivity (molecular and subgrid eddy). These 
are provided by the solution of Eqs. (6) -(12) by a finite 
difference LES (as described below) on a fixed grid and then 
interpolated to the particle location. 

The compositional values are subject to change due to 
subgrid mixing and chemical reaction. Equation (38) may be 
integrated numerically to simulate these effects simulta- 
neously. Alternately, this equation is treated in a split man- 
ner. This provides an analytical expression for the subgrid 
mixing. Subsequently, the influence of chemical reaction is 
determined by evaluating the fine grained reaction rates ui n a 
and modifying the composition of the elements. The imple- 
mentation of the SGS mixing and chemical reaction requires 
the filtered mean values of the scalars. These mean values 


(and other higher moments of the FDF) at a given point are 
estimated by consideration of particles within some volume 
centered at the point of interest Effectively, this finite vol- 
ume constitutes an “ensemble domain” characterized by the 
length scale A £ (not to be confused with A c ) in which the 
FDF is represented discretely by stochastic panicles. This is 
necessary as, with probability one, no particles will coincide 
with the point as considered. 56 Here, a box of size A £ is used 
to construct the ensemble mean values at the finite difference 
nodes. These values art then interpolated to the particle po- 
sitions. Since the mixing model only requires the input of the 
filtered scalar value, and not its derivative, this volume av- 
eraging procedure is sufficient However, from the numerical 
standpoint determination of the size of the ensemble domain 
is an important issue. Ideally, it is desired to obtain the sta- 
tistics from the Monte Carlo solution when the size of 
sample domain is infinitely small (i.e., A £ — ►O) and the num- 
ber of particles within this domain is infinitely large. With a 
finite number of particles, if A £ is small there may not be 
enough particles to construct the statistics. A larger ensemble 
domain decreases die statistical error, but may increase the 
dispersion error which manifests itself in “artificially dif- 
fused” statistical results. This compromise between the sta- 
tistical accuracy and dispersive accuracy as pertaining to La- 
grangian Monte Carlo schemes implies that the optimum 
magnitude of A £ cannot, in general, be specified a priori. 40 
This does not diminis h the capability of the procedure, but 
exemplifies the importance of the parameters which govern 
the statistics. 

The LES of the hydrodynamic variables, which also de- 
termines the subgrid viscosity and scalar diffusion coeffi- 
cients, is conducted with the “compact parameter” finite 
difference scheme of Carpenter. 77 This is a variant of the 
McCormack 73 scheme in which a fourth order compact dif- 
ferences are used to approximate the spatial derivatives, and 
a second order symmetric predictor-corrector sequence is 
employed for time discretization. The computational scheme 
is based on a hyperbolic solver which considers a fully com- 
pressible flow. Here, the simulations are conducted with a 
low Mach number (A/~ 0.3) to minimize compressibility ef- 
fects. The procedure involved in the finite difference discreti- 
zation is independent of the Monte Carlo solver, thus alter- 
native differencing schemes can be used if desired. All the 
finite difference operations are conducted on fixed and 
equally sized grid points. The transfer of information from 
these points to the locations of the Monte Carlo particles is 
conducted via interpolation. Both fourth-order and second- 
order (bilinear) interpolation schemes were considered, but 
no significant differences in SGS statistics were observed. 
The results presented in the next section are based on simu- 
lations with fourth- and second-order interpolations in two- 
dimensional (2D) and 3D flows, respectively. 

VI. RESULTS 
A. Rows simulated 

To demonstrate the effectiveness of the FDF method, in 
this section simulation results are presented of a temporally 
developing mixing layer and a spatially developing planar 
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jet. Both non-reacting and reacting flows are considered. In 
the latter, a simple reaction of the type is con- 

sidered in which the reaction is assumed to be constant rate 
and non-heat releasing (isothermal flow). Therefore, 

= o>#= — KAB , where K is a constant and A f B denote the 
mass fractions of species JB, respectively. The species 

Jff, are assumed thermodynamically identical and the 
fluid is assumed to be calorically perfect. Both 2D and 3D 
simulations are conducted of the temporal mixing layer. The 
jet simulations are 2D. 

The temporal mixing layer consists of two coflowing 
streams traveling in opposite directions with the same 
speed. 74 ” 77 The reactants and JB are introduced into the 
top and the bottom streams, respectively. In the planar jet, 
the reactant is injected with a high velocity from a jet of 
width D into a coflowing stream with a lower velocity car- 
rying reactant J?. 76,78 Both these flows are dominated by 
large scale coherent structures. The formation of these struc- 
tures are expedited by imposing low amplitude perturbations. 
In the figures presented below, x,y correspond to the stream- 
wise and cross-stream directions, respectively. In 3D, z de- 
notes the spanwise direction. In the temporal mixing layer, 
the length in the stream wise direction is chosen to be twice 
the wavelength of the most unstable mode to allow for the 
rollup of two large vortices and one (subsequent) pairing of 
these vortices. In 3D, the lengths in the streamwise and the 
cross-stream directions are the same as those in 2D. The 
length in the spanwise direction is 60% of that in the stream- 
wise direction. The forcing is imposed in such a way to 
provide significant 3D transport 79 ’ 80 The initial values of the 
mass fractions of reactants and JB at each of the spanwise 
points in 3D are identical to those in 2D. The size of the 
domain in the jet flow is 0^x^l4£>, — 3.5D^y^3.5D. 
The velocity ratio of the coflowing stream to that of the inlet 
jet is kept fixed at 0.5. 

Both flows are simulated via both DNS and LES. The 
procedure in DNS is exclusively based on the finite- 
difference solution of Eqs. (l)-(4) in which there are suffi- 
cient grid points to resolve the flow without a need for sub- 
grid closures. The procedure in LES is based on the Monte 
Carlo solution of the modeled FDF transport equation (Eq. 
(32)) for the scalars augmented by the finite difference solu- 
tion of the modeled equations of the filtered hydrodynamic 
variables (Eqs. (6)-(7)). In the presentation below, these re- 
sults are identified by the abbreviation FDF. In addition, an- 
other LES is conducted in which the modeled transport equa- 
tions for the filtered scalar and the generalized subgrid 
variance are simulated with the finite difference scheme. In 
these simulations, the hydrodynamic solver and the models 
for the subgrid stresses and mass fluxes are identical to those 
employed in FDF, but the effects of SGS fluctuations in the 
filtered reaction rate are ignored. Effectively, Eqs. (33)— (34) 
are solved with the assumption (w a (4>)) i: = <o a ((^) L ). The 
results based on this procedure are referred to as LES-FD. 
(The approximation was also con- 

sidered but did not show an improvement over LES-FD.) 

In both FDF and LES-FD simulations, the subgrid 
stresses are modeled via the Smagorinsky closure (Eqs. (9)- 
(10)) and the MKEV model (Eq. (11)). The subgrid mass 


fluxes are modeled via Eq. (12). No attempt is made here to 
determine the magnitudes of the constants appearing in these 
models in a dynamic manner. 14 However, several different 
values are considered for C, and C k . The values which give 
the best overall agreement with DNS in non-reacting flows 
are C x =0.014, 0.01 and C*=0.02, 0.013, in 2D, 3D, respec- 
tively. These values are subsequently used in FDF and 
LES-FD of scalar quantities in reacting flows. This param- 
eterization is justified since the LES of the hydrodynamic 
field is not the subject of our FDF closure. The other param- 
eters Sc= 1, Sc, =0.7 are kept fixed. In the MKEV model, 
the ratio of the filter size at the secondary level to that at the 
grid level is Ac»/A(;=3. In the implementation of the 
MKEV in the shear flows as considered, the magnitude of 
the reference velocity is set to zero in the cross-stream 
direction and to the average of the high and low speed 
streams in the streamwise direction. The subgrid mixing 
model requires the input of the constant in the mixing 
frequency which also determines the SGS variances. As will 
be shown below Cq *»3 is suggested, but the influence of 
this parameter is also studied by considering other val- 
ues. 

The flow variables are normalized with respect to refer- 
ence quantities denoted by the subscript r. In the temporal 
mixing layer the reference quantities are the freestrcam val- 
ues and the length L r is defined such that (S v0 fL r ) = 2.83, 
where S v0 is the initial vorticity thickness. In the planar jet, 
L r =D and the reference quantities are those at the high 
speed jet stream. The reference quantities define the Rey- 
nolds number Re= (f/^/v). For the temporal mixing layer, 
the Reynolds number based on the total velocity difference 
across the layer ( AU=2U r ) is given by Re^ 0 =5.66 Re. The 
reaction rate is parameterized by the Damkohler number 
Da =Kf(U r fL r ). The non-dimensional time is given by t* 
— (U r tfL T ). In the presentations below, the asterisk is 
dropped. 

B. Numerical specifications 

The magnitude of the flow parameters considered in the 
simulations are dictated by the resolution which can be af- 
forded by DNS. The primary parameters are the flow Rey- 
nolds number (Re), the Damkohler number (Da) and the mo- 
lecular Schmidt number. All finite difference simulations (in 
both DNS and LES) are conducted on equally-spaced, square 
(Ajc= Ay = Az (for 3D) = A) grids. Since the size of the 
computational domain is fixed, the number (and the size) of 
the grids depends on type of simulation being conducted. 
The highest resolution in DNS of the 2D temporal mixing 
layer consists of 433X577 grid points which allows reliable 
DNS with Re =500 and Da= 2 (based on the velocity differ- 
ence and the vorticity thickness at the initial time). The DNS 
of the 3D temporal mixing layer is conducted on 217X289 
X 133 grid points with Re =400, Da= 1. The resolution in 
DNS of the planar jet consists of 721X361 grid points and 
allows accurate simulations with Re= 12 000 and Da=2 
(based on the centerline velocity at the inlet and the jet 
width). 

The FDF and LES-FD are conducted on grids coarser 
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than those in DNS. Unless otherwise specified, the LES reso- 
lutions in the mixing layer consist of 37X49 grid points in 
2D, and 55X73X34 grid points in 3D. For the planar jet, a 
resolution of 101X51 is used for nonreacting flow simula- 
tions with Re=5000, while a 181X91 grid is utilized for 
reactive flow simulations with Da =2 and Re= 12 000. A 
top-hat filter function 47 of the form below is used 


G(*'-x)=n 

i-1 


(40) 


<7(x'-x,) = 




FIG. 1. 2D mixing layer sunuUnon results: Contours of the filtered con- 
served scalar, (a) FDF and (b) LES-FD. 


in which N D denotes the number of dimensions, and A c 
= 2A. No attempt is made to investigate the sensidvity of 
the results to the filter function 48 or the size of the filter. 81 

In FDF, the Monte Carlo particles are distributed at t 
= 0 throughout the domain. In the temporal mixing layer, the 
particles are distributed evenly throughout the whole compu- 
tational region. In the FDF of the jet, the particles are sup- 
plied initially in the inlet region — 1.75Z)«Sy^ 1.75D. In all 
the simulations, the particle density is monitored at all times 
to ensure an approximately uniform distribution throughout 
the mixing regions. In the temporal mixing layer, due to flow 
periodicity in the stream wise direction, if the particle leaves 
the domain at the right or the left boundary, new particles are 
introduced at the other boundary with the same composi- 
tional values. A similar procedure is employed in the span- 
wise direction in 3D simulations. Due to mirror symmetry at 
the upper and lower boundaries, particles that exit the top or 
bottom boundaries return to the domain at the opposite 
boundary with the mass fractions values associated with 
and JB interchanged. In the spatial jet, new particles are in- 
troduced at the inlet at a rate proportional to the local flow 
velocity and with a compositional makeup dependent on the 
y coordinate. The density of the Monte Carlo particles is 
determined by the initial number of particles per grid cell 
(NPG) of dimension AxA (xA in 3D). The magnitude of 
NPG is varied to assess its affect on statistical convergence 
of the Monte Carlo results. This assessment is demonstrated 
in 2D simulations of the temporal mixing layer. The simula- 
tions of 3D temporal layer and the spatial jet are based on 
NPG =20. The size of the “ensemble domain” in the FDF 
simulations is also varied to assess its influence on the sta- 
tistical convergence. The following sizes are considered: 
A f =2A,A,A/2. The number of samples used to construct 
the FDF is thus controlled by the values of NPG and A E . 

An additional parameter which influences the numerical 
accuracy is the magnitude of the incremental time step. The 
stability criterion for the finite difference scheme requires 72 
CFL=£ 1A/3 and is more stringent than the criterion for the 
Fourier number. The effect of the time increment on the 
accuracy of the Euler- Maruyamma scheme is also consid- 
ered. This is assessed by considering several At values (CFL 


numbers). In the context of Ito calculus, 82 ^ 3 this issue is 
considered by analysis of moments of the FDF up to the 
second order. 

The simulated results are analyzed both “instanta- 
neously” and “statistically.” In the former, the instanta- 
neous contours (snap-shots) of the scalar values are consid- 
ered. In the latter, the “Reynolds-averaged” statistics 
constructed from the instantaneous data are considered. In 
the spatially developing jet flow this averaging procedure is 
conducted via sampling in time. In the temporal mixing 
layer, the flow is homogeneous in x (and z in 3D); thus the 
statistics are generated by constructing the ensemble from all 
the grid points in the stream wise (and span wise) directions. 
These statistics are y — t dependent All Reynolds averaged 
results are denoted by an overbar. 

C. Consistency of FDF and convergence of the Monte 
Carlo procedure 

The objective in the results presented in this subsection 
is to demonstrate the consistency of the FDF formulation and 
the convergence of the Monte Carlo simulations. For this 
purpose, the LES results via FDF and LES-FD are compared 
against each other in shear flows under different conditions. 
In non-reacting flows, any deviations between the FDF and 
LES-FD results are attributed to the differences in the nu- 
merical procedures. Since the accuracy of the finite differ- 
ence procedure is well-established, this comparative analysis 
provides a good means of assessing the performance of the 
Monte Carlo solution of the FDF. Unless specified other- 
wise, the Smagorinsky model is used to evaluate the eddy 
viscosity in the simulations considered in this subsection. 

In Fig. 1, results are presented of the LES of the non- 
reacting temporally developing mixing layer. Shown in the 
figure are the contour plots of (A) L via (a) FDF and (b) 
LES-FD, with A = 0, 1 on the bottom and top streams, re- 
spectively. These contours correspond to results at a time 
when the flow has experienced the pairing of two neighbor- 
ing vortices. This figure provides a simple visual demonstra- 
tion of the consistency of the FDF as the results via the 
particle method are in agreement with those obtained by 
LES-FD. In fact, the Monte Carlo results are somewhat more 
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attractive due the Lagrangian nature of the solution proce- 
dure. While the LES-FD results display slight over- and 
under-shoots, there are no such errors in the Monte Carlo 
scheme. 

A more rigorous means of assessing the FDF results is 
via consideration of the Reynolds averaged results. Figures 2 
and 3 show such results in the non-reacting temporal mixing 
layer in which the sensitivity of the FDF predictions to sev- 
eral parameters is assessed. Figure 2(a) sh ows the compari- 
son of FDF and LES-FD results for (A) L for several values 
of A £ . It is shown that the first filtered moment of the FDF 
agrees very well with that obtained by LES-FD, even for 
large A £ values. The differences between the FDF and 
LES-FD results are more apparent in Figs. 2(b,c,d) where the 
cross-stream variations of cr A are shown for several values of 
A £ and and for different LES grid resolutions. As ex- 
pected, Figs. 2(b,c) show that with increasing C*n , the mag- 
nitude of the variance decreases. These figures also indicate 
that the difference between FDF and LES-FD predictions 
diminish as A £ decreases. This is also corroborated in Fig. 
2(d) in which the both FDF and LES-FD are conducted on 
61 X 81 grid points. At all A £ values, the agreement between 
FDF and LES-FD is better than those shown in Fig. 2(b) 
with a lower finite difference resolution. The consistency of 
the FDF and LES-FD results does not mean that the magni- 
tude of Cn can be specified. Hereinafter Cn=3 is adopted 
since it provides the best overall match with DNS data as 
shown in the next subsection. 

The other parameters which influence the accuracy of 
the Monte Carlo results are the number of Monte Carlo par- 
ticles per grid cell (NPG) and the magnitude of the incre- 
mental time step. Figure 3(a) shows that a A values do not 
vary significantly for NPG>50. In fact in other cases even 
smaller NPG values can be used as will be shown. Figure 
3(b) verifies the insensitivity of statistics to A t as long as the 
stability criterion is satisfied (CFL<l/vG). Hereinafter, 
CFL=0.5 is used. 

The sensitivity of the results to NPG and A E in the FDF 
simulations of a reacting temporal mixing layer with Da= 2 
is studied in Fig. 4. In these simulations, the MKEV model is 
adopted to evaluate the subgrid viscosity because it performs 
somewhat better than the Smagorinsky model for LES of 
reactive flows (as assessed by DNS data in the next subsec- 
tion). Shown in the figure are the Rey nolds averaged values 
of the filtered product mass fraction ((P) £ ) at a fixed time 
(Fig. 4(a)) and the integrated total product (8 P (t) 
= f{P)i(yJ)dy). The convergence of Monte Carlo solution 
and the independence to NPG and A E arc demonstrated by 
these results (at least for this first moment). Moreover, it is 
shown that while the mean value of the scalar as used in the 
mixing model for a given particle should be evaluated at the 
particle location, the mean value at the nearest finite differ- 
ence grid point could also be substituted. This eliminates the 
need for interpolating the mean scalar field to the particle 
locations. 

The consistency and the convergence of the Monte Carlo 
simulation of the FDF in the nonreacting jet flow arc sum- 
marized in Figs. 5-6 in which the time averaged (Reynolds) 






FIG. 2. 2D muring layer simulation results: (a) Conserved scalar distribution 
vs. the cross-stream coo rdinate Generalized variance vs. cross-stream coor- 
dinate at (b) C n =l and (c) C n « 3. (d) Same as (b) but with increased 
resolution. 
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FIG. 3. 2D mixing layer simulation results: (a) Cross-stream variation of the 
generalized variance for various NPG. (b) Generalized variance vs. cross- 
stream coordinate at various CFL numbers. 




FIG. 4. 2D mixing layer simulation results: (a) Cross-stream variation of the 
product mass fraction, (b) Total product vs. trine. The long-dashed line 
represent the case where the mean value of the scalar in the mixing model 
for a particle is set to be equal to the value at the nearest finite difference 
grid point. 



FIG. 5. 2D planar jet results: Conserved scalar distribution vs. cross-stream 
coordinate at x~5D r 9D. 


stat is tics for the first and second subgrid moments of A are 
presented. Similar to the temporal mixing layer results. Fig. 
5 shows the accuracy of the Monte Carlo solver and the 
insensitivity of results to A £ for the first moment of the FDF. 
Si milar ly, for the scalar variance, the agreement between the 
FDF and LES-FD results diminishes as the size of A £ is 
decreased. At x=5D, the FDF results with A £ = A are very 
close to those via LES-FD. With the same A £ values the 
agreement is not as good at x — 9D and lower values of A £ 
are needed to achieve a better agreement for the subgrid 
variance. However, as will be shown below, with this reso- 




F1G. 6. 2D planar jet results: Generalized variance vs. cross-stream coordi- 
nate at (a) x=5 D and (b) x~9D. 
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RG. 7. 3D mixing layer simulation results: Scatter plot of the filtered values 
of a conserved scalar as obtained by FDF vs. those via LES-FD. 

lution the mean filtered values of reacting scalars are pre- 
dicted reasonably well. 

The consistency of the FDF simulation in 3D is demon- 
strated in Fig. 7 in which the scatter plot is shown of the 
instantaneous filtered A values as obtained by FDF vs. those 
via LES-FD. The hydrodynamic LES is based on MKEV in 
both simulations. The correlation coefficient between the 
data obtained by the two simulations is 0.99. This excellent 
correlation is also reflected in the cross stream profiles of the 
Reynolds-averaged filtered quantities in Fig. 8. 

D. DNS validations of the FDF 

The objective in this section is to assess the overall per- 
formance of the FDF and to appraise the validity of the sub- 
models employed in the FDF transport equation. For this 
objective, the FDF results are compared against DNS of the 
same flow configuration with the same magnitudes of Re and 
Da. For a meaningful comparison, the DNS data are filtered 
and the results on the coarse grids are compared with those 
on the corresponding coarse grids in the FDF simulations. To 
illustrate the capability of the FDF, the results are also com- 
pared with LES-FD in which the effects of SGS fluctuations 
on the filtered reaction rate are ignored. 




(b) 


RG. 9. 2D mixing layer simulation results: Effect of grid resolution on 
temporal evolution of the (a) vortidty thickness and (b) total product. 


First the resolution requirement for DNS is determined. 
This is demonstrated here for the 2D mixing layer. A similar 
procedure is followed for the other flow configurations. In 
Fig. 9 results are presented of the temporal evolution of the 
vorticity thickness (£„) and the total product ( S P ) in a react- 
ing layer with Re =500, Da= 2 at several resolutions. It is 
shown that the results generated via 289x385 are almost 
identical to those on 433X577 grid points. Analysis of other 
statistical results (not shown) show a similar behavior. Here- 
inafter 433 X 577 grid points are used in all DNS of the 2D 
mixing layer. The resolution employed in I-F.S (both FDF 



RG. 8. 3D mixing layer simulation results: Cross-stream variations of the RG. 10. 2D mixing layer simulation results: The integrated Reynolds aver- 
mean value of the filtered mass fraction of a conserved scalar. aged values of the filtered scalar's va riance vs. time. 
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FIG. 11. 2D mixing layer simulation results: Total SGS unnrixedness and 
Reynolds subpan vs. cross-stream coordinate. 

and LES-FD) is coarser consisting of 37X49 grid points. 
The results in Fig. 9 indicate the inaccuracy of “DNS” at 
this resolution. 

To determine the magnitude of C n , in Fig. 10 the inte- 
grated Reynolds averaged values of the SGS variance 
(S^A(y^)dy) of a nonreacting scalar as predicted by FDF 
are compared with those via DNS. This comparison shows 
that C n — 3 yields a reasonable agreement between the pre- 
diction and DNS results. Thus, this value is used in absence 
of a better model of the subgrid mixing frequency. 

To demonstrate the difficulty of modeling the SGS scalar 
fluctuations in reacting flows, the Reynolds averaged profiles 
for the “SGS unmixedness” (t ab = (AB) l — (A) l (B) l ) and 
its “Reynolds” subpart 84 * 85 9 )i(B 9 ) l as 

obtained directly from DNS data are shown in Fig. 11. These 
results show the importance (non-zero values) of these cor- 
relations. They also show that R AB is a fraction of r ^ sug- 
gesting that modeling of t ab in LES is more complex than 
that in Reynolds procedures. 

In Fig. 12, the FDF predictions of the total product are 
compared with DNS results. The Smagorinsky model is em- 
ployed in FDF with several values of the parameter C s . Ob- 
viously for a constant C s value, the agreement between DNS 
and FDF is not very satisfactory. The subgrid viscosity based 
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FIG. 12. 2D mixing layer simulation results: Total product variation with 
time. The Smagonnsky model is used to represent the eddy viscosity for the 
FDF simulations. 



FIG. 13. 2D mixing layer simulation results: Voracity thickness vs. time. * 

on the Smagonnsky closure affects both the resolved hydro- 
dynamic field and the subgrid scalar mixing process. It is 
known that the Smagonnsky closure sometimes generates 
excessive damping on the resolved scales in transitional 
regions. 49 Here, an attempt is made to rectify the situation, 
albeit in a very ad hoc manner. In the temporal mixing later, 
C t should be initially zero to reflect the fact that the flow is 
“laminar.” Then its value should increase in time as the 
flow becomes more “turbulent” The FDF results in Fig. 12 
with Cj^t agree more favorably with DNS. This is partly 
due to better predictions of the hydrodynamic field (Fig. 13) 
but also due to more accurate representation of the subgrid 
mixing frequency. This better agreement is not sufficient to 
suggest a new model for C, ; rather it is to demonstrate the 
importance of the subgrid diffusion in affecting the FDF di- 
rectly (through the subgrid mixing) and indirectly (through 
the input of the hydrodynamic parameters). 

In order to better predict the subgrid viscosity, the 
MKEV model (Eq. (11)) is adopted. In Fig. 13 it is shown 
that the vorticity thickness predicted by the MKEV model 
compares with DNS data better than that via the Smagorin- 
sky model. The improved prediction of the eddy viscosity 
also improves the FDF predicted product formation as shown 
in Fig. 14 for several values of the Damkohler number. Due 
to the demonstrated superiority, the MKEV closure is uti- 
lized in all subsequent simulations unless otherwise noted. 



FIG. 14. 2D mixing layer simulation results: Temporal evolution of the total 
product. 





FIG. 16. 2D mixing layer simulation results: Scatter plots of instantaneous 
value of the conserved scalar vs. the mean value. Data taken from (a) DNS, 
(b) FDF throughout the computational domain. 


FIG. 18. 2D planar jet simulation results: Instantaneous reaction rate as 
determined by (a) DNS. (b) FDF, (c) LES-FD. 
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FIG. 19. 2D planar jet simulation results: Cross-stream variation of the 
mean product mass fraction at (a) x = 5D and (b) x*9D. 


To demonstrate the importance of the SGS scalar fluc- 
tuations, the results of FDF and LES-FD are compared with 
DNS results in Fig. 15. This figure shows that the neglect of 
SGS unmixedness results in significant overpredictions of 
the product mass fraction. This behavior is observed at all 
times and all values of the Damkohler number (Fig. 15(b)) 
and is consistent with that in Reynolds averaging. 18 More- 
over, Fig. 15(b) shows that as the magnitude of the 
Damkohler number increases, the neglect of the SGS unmix- 
edness in LES-FD results in progressively higher deviation 
of product formation relative to DNS. This is significant 
since the Da values in practical reacting systems can be quite 
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FIG. 21. 3D mixing layer simulation results: Cross-stream variation of the 
product distribution. 


large. Therefore it is expected that the effects of the SGS 
unmixedness are very pronounced in such applications. To 
verify that the enhanced product formation in LES-FD is not 
associated with the numerical discretization errors, an addi- 
tional FDF is conducted in which the filtered reaction rate is 
“incorrectly” calculated in termsof the filtered values of the 
reactants' mass fractions. The results based on this model are 
identified by FDF* in Fig. 15(a) and consistent with LES-FD 
results, overpredict the rate of reactants' conversion. 

It is useful to compare the DNS results for “fine grid” 
scalar values with the “fine-grained” values associated with 
the Monte Carlo particles. The “scatter” plots of the instan- 
taneous fine grid values of A vs. its filtered value {A) L as 
obtained by DNS are presented in Fig. 16(a) and the scatter 
plot of fine grained A values vs. (A) L is shown in Fig. 16(b). 
These results are associated with a non-reacting temporal 
mixin g layer and are taken at a fixed time. The points in Fig. 
16(a) correspond to the values at all the grid points employed 
in DNS within the computation domain. The points in Fig. 
16(b) correspond to all Monte Carlo panicles occupying the 
same domain. It is shown that the “density” of scatter is 
similar in the two plots indicating a qualitative agreement 
between FDF and DNS. However, the scatter in FDF is ex- 
pectedly somewhat greater but not with a significant density. 



t 


FIG. 20. 2D planar jet simulation results: Total product vs. the downstream 
coordinate. 


FIG. 22. 3D mixing layer simulation results: Temporal evolution of the total 
product. 
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TABLE I. Total computational times for the 2D reacting mixing layer simu- 
lations. 


Simulation 

Grid resolution 

NPG 

Normalized CPU time* 

Figure 

DNS 

433X577 

— 

285.45 

14, 15(b) 

FDF 

37X49 

40 

8.45 

14 

LES-FD 

37X49 

— 

1 

15(b) 


•Unit corresponds to 11 s on a Cray-C90. 


The effectiveness of the FDF to predict the slightly more 
complex jet flow is summarized in Figs. 17-20. Figure 17 
shows the instantaneous contours of the normalized SGS un- 
mixedness as obtained by filtered DNS and FDF. Note that 
this term is assumed to be identically zero in LES-FD. The 
SGS unmixedness is negative throughout the reaction zone, 
thus its effect is manifested in a decrease of the filtered re- 
action rate. This is readily observed in Fig. 18, where the 
contour plots of the reaction rates are displayed for the fil- 
tered DNS, FDF and LES-FD approaches. While the peak 
values in the DNS are slightly higher than those observed in 
the FDF simulations, the reaction zone predicted by the FDF 
simulation is slightly thicker (due to the finite size of the 
ensemble domain) therefore yielding a comparable amount 
of converted products. In contrast* the filtered reaction rates 
obtained by the finite difference LES procedure in which the 
SGS unmixedness is neglected are significantly higher. This 
is reflected in Fig. 19, where the cross-stream variation of the 
time-averaged filtered values of the product mass fraction are 
presented at two downstream locations and in Fig. 20, where 
the strea mwise variation of the integrated total product 
(S P (x) = J(P)L(x>y)dy) is shown. Two additional points 
are intended by presentations of Figs. 19 and 20. First, the 
FDF results are compatible with those of DNS at all down- 
stream coordinates. Therefore, there is no “secular” behav- 
ior associated with possible modeling err o r s in the FDF. Sec- 
ond, the differences between the FDF and DNS in predicting 
the subgrid scalar variances at large xfD values as observed 
in the variance results in Fig. 6 do not seem to yield signifi- 
cant differences in the amount of product formation as pre- 
dicted by the FDF. In all the cases the neglect of the SGS 
fluctuations, as done in LES-FD, results in significant over- 
predictions of the filtered reactant conversion rate. It is ex- 
pected that these overpredictions would become even more 
significant at higher Damkohler and Reynolds numbers. 

The major conclusions drawn from the 2D results are 
confirmed in 3D simulations. The cross-stream variation of 
the filtered mean products and the temporal variation of the 
total product in the 3D mixing layer are shown in Figs. 21 
and 22. The performances of the Smagorinsky and MKEV 


TABLE n. Total computational times for the reacting jet simulations. 


Simulation 

Grid resolution 

NPG 

Normalized CPU tune* 

Figure 

DNS 

721X361 

— 

52.12 

18(a) 

FDF 

181X91 

20 

12.56 

18(b) 

LES-FD 

181X91 

— 

1 

18(c) 


•Unit corresponds to 809 s on a Cray-C90. 


TABLE IR Total co m p uta tional times for the 3D reacting mixing layer 
simulations. 


Simulation 

Grid resolution 

NPG 

Normalized CPU time* 

Figure 

DNS 

217X 289X133 

_ 

182.71 

21.22 

FDF 

55X73X34 

20 

7.64 

21,22 

LES-FD 

55 X 73 X 34 

- 

1 

21.22 


•Unit correspond to 655 s on a Cr*y-C90. 


closures in predicting the hydrodynamic field are similar to 
those in 2D. With either closures, the amount of products 
predicted by LES-FD is higher than those obtained by FDF 
and DNS. The FDF results are again in a good agreement 
with DNS flat* This agreement also indicates that the mix- 
ing model with Cfl = 3 works well in 3D simulations; no 
attempt was m aflg to find the optimize value of this constant. 
Future applications to other flow configurations would deter- 
mine the generality of the model. 


E. Comparison of computational requirements 

The total computational times associated with some of 
the simulations are shown in Tables I— ID. The cases consid- 
ered in this table are those which give reasonably accurate 
predictions of the first FDF moments of the reacting scalar 
field. Expectedly, the overhead associated with the FDF 
simulation is somewhat extensive as compared to LES-FD; 
nevertheless the FDFs computational requirement is signifi- 
cantly less that of DNS. While this overhead was tolerated in 
present simulations, there are several means of reducing it 
for future applications. A detailed examination of the indi- 
vidual routines utilised in the FDF simulations indicates that 
the most demanding computation is associated with the par- 
ticle interpolation procedure. The fourth order interpolation 
routine consumes about 51.3% of the total CPU time. The 
bilinear scheme reduces the computational time by 36%. If 
interpolation can be totally disregarded, i.e., using the results 
at the nearest finite difference grid point as shown in Fig. 4, 
the CPU time can be decreased by 50%. In addition, the 
Lagrangian procedure would benefit from the utilization of 
parallel architecture, since a significant portion of the time is 
devoted to computations in large loops dimensioned by the 
total number of Monte Carlo particles. This has been dis- 
cussed for use in PDF (Ref. 86) and its utilization is recom- 
mended for FDF. 

In comparing the computational requirements of FDF 
with those of DNS, it is important to note that this compari- 
son could be made only in flows for which DNS was pos- 
sible. The DNS times and the FDF times are as close as they 
are simply because the DNS had to be done at low Re, Da 
values. At higher values of these parameters, the difference 
could be much greater. This warrants further extensions and 
applications of FDF for more complex turbulent reacting 
flows for which DNS is not possible. 
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VII. CONCLUDING REMARKS 

It is demonstrated that the filtered density function 
(FDF) provides a powerful method for large eddy simulation 
(LES) of turbulent reacting flows. The method is based on 
the representation of the distribution of the unresolved fluc- 
tuations at the subgrid level. This is similar to the probability 
density function (PDF) methods in Reynolds averaging pro- 
cedures. Here, the FDF methodology is developed for treat- 
ment of scalar variables. Thus, similar to PDF methods it 
represents the effects of chemical reactions in a closed form. 

A modeled transport equation is developed for the FDF 
by adopting a closure strategy similar to that in PDF meth- 
ods. It is shown that the Lagrangian Monte Carlo scheme 
provides an effective means of solving the FDF transport 
equation. The scheme is exploited for I ES of two- and three- 
dimensional shear flows under both nonreacting and reacting 
conditions. The simulated results are compared with those 
based on conventional LES methods in which the effects of 
subgrid scalar fluctuations are ignored ( LES -FD), and those 
via direct numerical simulation (DNS) of flows with identi- 
cal values of the physical parameters. The convergence of 
the Monte Carlo numerical results and the consistency of the 
FDF formulation are demonstrated by comparisons with the 
Eulerian results of LES-FD of non-reacting flows. The supe- 
riority of the FDF over LES-FD is demonstrated by detailed 
comparative assessments with DNS results of reacting shear 
flows. It shown that the subgrid scale scalar fluctuations have 
a very significant influence on the filtered reaction rate; the 
neglect of these fluctuations results in overpredictions of the 
filtered reactant conversion rate. 

Although the present methodology is developed for iso- 
thermal, constant density, reacting flows with a simple kinet- 
ics scheme, the extension to variable density flows, with exo- 
thermic reactions imposes no serious mathematical 
difficulties. With such an extension, it is conceivable that 
LES of reactive flows with realistic chemical kinetics may be 
conducted for engineering applications in the near future, if 
the computational overhead associated with the FDF can be 
tolerated. In this regard, the scalar FDF methodology is at- 
tractive in that the present Monte Carlo solver can be used 
directly in available CFD codes. Similar to PDF methods, the 
closure problems associated with the FDF are the correla- 
tions involving the velocity field (such as SGS stresses and 
mass fluxes). This may be overcome by considering the joint 
velocity-scalar FDF similar to that in PDF. 87 

The computational requirement for FDF is more than 
that for LES-FD and less than that for DNS. The range of 
flow parameters (such as the Reynolds and the Damkohler 
numbers) that can be considered by FDF is significantly 
larger than can be treated by DNS, and the results are more 
accurate that those by LES-FD. This comparison of compu- 
tational requirements could be made here only in flows for 
which DNS was possible, i.e., low Da, Re values. At higher 
values of these parameters, the computational cost associated 
with DNS would be exceedingly higher than that of FDF. 
Thus for practical flows for which DNS is currently impos- 
sible, FDF would be a good alternative. Several means of 
reducing the FDF s computational requirements are recom- 


bia 

mended. These could be useful in future applications in com- 
plex flows. The FDF methodology will benefit from ongoing 
and future improvements in PDF schemes from both model- 
ing and computational standpoints 56 
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Abstract 

A methodology termed the “filtered mass density function” (FMDF) is developed 
and implemented for large eddy simulation (LES) of variable density chemically react- 
ing turbulent flows at low Mach numbers. This methodology is based on the extension 
of the “filtered density function” (FDF) scheme recently proposed by Colucci et al. 
(1998) for LES of constant density reacting flows. The FMDF represents the joint 
probability density function of the subgrid scale (SGS) scalar quantities and is ob- 
tained by solution of its modeled transport equation. In this equation, the effect of 
chemical reactions appears in a closed form and the influences of SGS mixing and 
convection are modeled. The stochastic differential equations (SDEs) which yield sta- 
tistically equivalent results to that of the FMDF transport equation are derived and 
are solved via a Lagrangian Monte Carlo scheme. The consistency, convergence, and 
accuracy of the FMDF and the Monte Carlo solution of its equivalent SDEs are as- 
sessed. In non-reacting flows, it is shown that the filtered results via the FMDF agree 
well with those obtained by the “conventional” LES in which the finite difference solu- 
tion of the transport equations of these filtered quantities are obtained. The advantage 
of the FMDF is demonstrated in LES of reacting shear flows with nonpremixed reac- 
tants. The FMDF results are appraised by comparisons with data generated by direct 
numerical simulation (DNS) and with experimental measurements. In the absence of 
a closure for the SGS scalar correlations, the results based on the conventional LES 
are significantly different from those obtained by DNS. The FMDF results show a 
closer agreement with DNS. These results also agree favorably with laboratory data of 
exothermic reacting turbulent shear flows, and portray several of the features observed 
experimentally. 
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1 Introduction 


Within the past decade, large eddy simulation (LES) of turbulent reacting flows has been 
the subject of widespread investigations (McMurtry et al . , 1992; Menon et al. , 1993; Gao 
and O’Brien, 1993; Madnia and Givi, 1993; Frankel et al., 1993; Cook and Riley, 1994; 
Fureby and Lofstrom, 1994; Moller et al., 1996; Branley and Jones, 1997; Cook et al., 1997a; 
Cook et al., 1997b; Jimenez et al., 1997; Mathey and Chollet, 1997; Colucci et al., 1998; 
DesJardin and Frankel, 1998; Jaberi and James, 1998; Reveillon and Vervisch, 1998); see 
Galperin and Orszag (1993); McMurtry et al. (1993); Libby and Williams (1994); Fox (1996); 
Vervisch and Poinsot (1988) for reviews. Amongst these, Colucci et al. (1998) recently 
developed a methodology, termed the “filtered density function” (FDF) based on an idea 
originally proposed by Pope (1990). The fundamental property of the FDF is to account for 
the effects of subgrid scale (SGS) scalar fluctuations in a probabilistic manner. Colucci et 
al. (1998) developed a transport equation for the FDF in constant density flows in which 
the effects of unresolved convection and subgrid mixing are modeled similarly to those in 
“conventional” LES, and Reynolds averaging procedures. This transport equation was solved 
numerically by a Lagrangian Monte Carlo procedure and the results were compared with 
those obtained by direct numerical simulation (DNS) and by a conventional finite difference 
LES in which the effects of SGS scalar fluctuations are ignored (LES-FD). It was shown that 
in non-reacting flows, the first two SGS moments of the FDF as obtained by the Monte Carlo 
solution are close to those obtained by LES-FD. The advantage of the FDF was demonstrated 
in reacting flows in which its results were shown to deviate significantly from those obtained 
by LES-FD but compare favorably with DNS data. 

The encouraging results generated by FDF warrant its extension and application to more 
complex flows. Further assessment of its predictive capability is also in order. The primary 
objective in this work is to extend the FDF methodology for treatment of variable density re- 
acting flows so that exothermic chemical reactions can be simulated. For that, we introduce 
the “filtered mass density function” (FMDF). With the definition of the FMDF, the math- 
ematical framework for its implementation in LES of reacting flows is established. A new 
computational scheme is also developed for the solution of the FMDF transport equation. 
The results obtained by FMDF are scrutinized by comparisons with DNS and laboratory 
data in several turbulent reacting flows with nonpremixed reactants. The FMDF deals only 
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with scalar quantities; the hydrodynamic field is obtained via conventional LES. Also, the 
formulation is based on the assumption of low Mach number. This allows consideration of 
exothermicity and variable density effects, but the method cannot be used for LES of very 
high speed flows (Drummond, 1991). 


2 Governing Equations 

In a compressible flow undergoing chemical reaction, the primary transport variables are the 
density p, the velocity vector it,, i = 1,2,3 along the x,- direction, the total specific enthalpy 
h, the pressure p, and the species mass fractions Y a {a = 1 , 2 ,..., N t ). The conservation 
equations governing these variables are the continuity, momentum, enthalpy (energy) and 
species mass fraction equations, along with an equation of state (Williams, 1985) 


dp dpui 

dt dxi 

= 0 

( 1 ) 

dpuj dpuiUj 

dt dxi 

dp drij 

dxj dx{ 

( 2 ) 

dp4> a dpui4> a dJf 0 

dt dxi ~~ dxi +p a ’ 

a = 1,2, . . . ,<r = N, + 1 

( 3 ) 

p = pR°T £ YJM q = pTZT 

( 4 ) 


0=1 


where t represents time, BP is the universal gas constant find Ai a denotes the molecular 
weight of species a. Equation (4) effectively defines the mixture gas constant 71. Equation 
(3) represents the transport of the species’ mass fractions and enthalpy in a common form 
with 

N t 

<f> Q =r Y a , a = 1,2,..., N s , <f> a = h = } " h a <f> a (5) 

0=1 

with 

K = hi + C C r „(T)dT' ( 6 ) 

JTq 

where T denotes the temperature, To is the reference temperature and and Cp a denote 
the enthalpy at To, and the specific heat of species a at constant pressure, respectively. At 
low Mach numbers and heat release rates, by neglecting the viscous dissipation and thermal 
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radiation the source terms in the enthalpy equation S a = fa can be assumed to 
be negligible. Thus, the chemical source terms ( S a = = [Vi, > 2 , . . . , Yn 3 , A]) are 

functions of the composition variables (4>). For a Newtonian fluid with zero bulk viscosity 
and Fickian diffusion, the viscous stress tensor r,j, mass and heat flux (J“, a = 1, 2, . . . , a) 


are given by 


Tij = H 


dui duj 
dxj **" dxi 


2du* 

Zdx k ij 




(V 

( 8 ) 


where fi is the dynamic viscosity and 7 = pT denotes the thermal and the mass molecular 
diffusivity coefficients. Both f* and 7 tire assumed constant and the Lewis number is assumed 
to be unity. In reactive flows, molecular processes are much more complicated than portrayed 
by Eq. (8). But since the molecular diffusion is typically less important than the SGS 
diffusion (to be defined below), this simple model is adopted with justifications and caveats 
given by Pope (1985); Bilger (1982). 


Large eddy simulation involves the use of the spatial filtering operation (Aldama, 1990) 


/ +00 

/(x',f)<y(x',x)dx' 

-OO 


( 9 ) 


where Q denotes the filter function of width Ac, {/(x,i))* represents the filtered value 
of the transport variable and f' = f— ( f) e denotes the fluctuations of / from the 

filtered value. In variable density flows it is convenient to consider the Favre filtered quantity 
(/( x > 1))l = (pf)e/(p)e and the fluctuation f" = /— {/)l- We consider spatially temporally 
invariant and localized filter functions, ^(x / ,x) = G(x' — x) with the properties (Aldama, 
1990), G(x) = G(-x), and G(x)dx = 1. Moreover, we only consider “positive” filter 
functions as defined by Vreman et al. (1994) for which all the moments x m G(x)dx exist 
for m > 0. The application of the filtering operation to the transport equations yields 


d(p)* ^(p)r(u,)L 

dt dx i 


( 10 ) 


djpUjujh , ^( p )/(^») l(^)l 

dt dxi 


d(p)> dim), dTj, 

dxj dxi dii 


( 11 ) 
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d(p)t(<f>a)L _ d(J°)e _ dM? 

dt dx i dxi dx{ 


+ {pS a )t, a = 1,2, ... ,cr (12) 


where = (p) t ((uiUj) L - (ui) L (uj) L ) and M? = (p)t({ui<f> a ) L - { Ui) L (<f> a ) L ) denote the 
subgrid stress and the subgrid mass flux, respectively. The filtered reaction source terms are 
denoted by ( pS a ) ( = ( p)t(S a ) L (a = 1 , 2 ,..., N t ). 


Modeling of Hydrodynamic SGS Quantities 


In LES of non-reacting flows the closure problem is associated with Tij and M° (Erlebacher 
et al., 1992; Salvetti and Banerjee, 1995). In reacting flows, an additional model is required 
for the filtered reaction rate ( S a )j_, . This is the subject of the probability formulation as 
described in the next section. For T \j, the variable density form of the model used in our 
previous work (Colucci et al., 1998) is considered: 


Tij = ~ 2C R (p)<A G £ 1/2 + IciWtSSii (13) 

where is the resolved strain rate tensor, S = |(«*)£,(w*)t — u i = 

Ui - Ui and U, is a reference velocity in the x, direction. The subscript i' denotes the filter 
at the secondary level of size Ac > Aq . This model is essentially a modified version of that 
proposed by Bardina et al. (1983), which utilize equal sizes for the grid and secondary filters. 
We refer to this as the modified kinetic energy viscosity (MKEV) closure. Accordingly, the 
subgrid eddy viscosity is expressed as v t = CrA g S^. A similar diffusivity model is used for 
the closure of the subgrid mass flux (Eidson, 1985) 


M Q = 


-It 


d(4>ah 

dxi 


(14) 


where 7 t = (p)eT t , I\ = v t /Sc t , and Sc t is the subgrid Schmidt number, assumed to be 
constant and equal to the subgrid Prandtl number. It must be emphasized here that these 
models are not used directly in the FMDF but the modeled FMDF transport equation is 
constructed to be consistent with them as discussed below. 
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3 The Filtered Mass Density Function (FMDF) 


Let <f>(x, t) denote the scalar array. We define the “filtered mass density function” (FMDF), 
denoted by Fl, as 


/ +oo 

/>(*', *)( tyb <£(x', *)] G{x! - x)dx', 

-oo 


(15) 


C [V>> ^(x, ()] = - <Hx, i)] - n < ('/'o - M x > <)] ( 16 ) 

a=l 

where 6 denotes the delta function and ip denotes the composition domain of the scalar array. 
The term ip(x, t)] is the “fine-grained” density (O’Brien, 1980; Pope, 1985), and Eq. 
(15) implies that the FMDF is the mass weighted spatially filtered value of the fine-grained 
density. The integral property of the FMDF is such that 

/ +oo r+oo 

F L {ip;x,t)dip = f p(x',t)G(x' - x)dx' = (p(x,t)) e . (17) 

-oo J — oo 


For further developments, the mass weighted conditional filtered mean of the variable Q(x, t) 
is defined as 


(Q(x,t)\ip) e 


M x V)Q( x '>0C[Vb<fr(x',*)] G{xf - x)dx' 
Fl{*P ; X, t) 


Equation (18) implies 


(18) 


(«') 

(«) 

(Hi) 


For Q(x, t) = c, ( Q(x , OIV’)/ = c 

For Q(x, t ) = Q{4 i>(x, t )), <<2(x, t)\ip) t = Q(ip) 

/ +oo 

(Q(x, t)\ip) t F L (ip; x, t)dip 

-OO 


(19) 

(20) 

( 21 ) 


where c is a constant, and Q(<p(x, t)) = Q(x,t) denotes the case where the variable Q can 
be completely described by the compositional variable <p(x , t ) = \<j> i, <f>z , . . . , <p a ]. From these 
properties, it follows that the filtered value of any function of the scalar vaxiables (such as 
p = p[<p(x,t)] and S a = 5 Q [<£(x,f)] ) is obtained by integration over the composition space. 
It is noted that the mass weighted conditional filtered mean reduces to the conditional filtered 
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mean (Colucci et al., 1998) when the density can be completely expressed in terms of the 
compositional variables. 

By applying the method developed by Lundgren (1969); Pope (1976); O’Brien (1980) to Eq. 
(3), a transport equation is obtained for the fine-grained density (Colucci et al., 1998). The 
transport equation for ; x, t) is obtained by multiplying the equation for the fine grained 
density by the filter function G(x' — x) and integrating over x' space. The final result after 
some algebraic manipulation is 

dF L (if>]*,t) d[(ui(x,t)\i/>) t F L (ijr,x,t)} __ d [7 1 dJ? \ 

— at — + = m^l, t{ *' ’ ' 

a&WFi.(^;x.0| 
dip a 

This is an exact transport equation for the FMDF. The last term on the right hand-side 
of this equation is due to chemical reaction and is in a closed form. The unclosed nature 
of SGS convection and mixing is indicated by the conditional filtered values. These terms 
are modeled in a manner consistent with Reynolds averaging and conventional LES in non- 
reacting flows. The convection term is decomposed via 

(23) 

where the second term on the right hand side denotes the influence of SGS convective flux. 
This term is modeled as 

(24) 

The advantage of the decomposition (Eq. (23)) and the subsequent model (Eq. (24)) is that 
they yield results similar to that in conventional LES (Germano, 1992; Salvetti and Banerjee, 
1995). The first Favre moments corresponding to Eqs. (23) and (24) are 

(Ui<t>a)L — (^i “1“ L fat) L($at} l\i (25) 

( 26 ) 


(ui\ip)tF L = {u^lFl + [(u.lV’)/ — (ui) L ]F L . 
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The term within brackets in Eq. (25) is the generalized scalar flux. This makes Eq. (26) 
identical to Eq. (14). The closure adopted for the SGS mixing is based on the linear mean 
square estimation (LMSE) model (O’Brien, 1980; Dopazo and O’Brien, 1976), also known 
as the IEM (interaction by exchange with the mean) (Borghi, 1988) 


d 

dx/} a 



s + JL _ (*„> 4 )*y , (27) 

OXi \ UXx J 


where D m (x, t) is the “frequency of mixing within the subgrid” which is not known a priori. 
This frequency is modeled as fi m = Cu{ 7 + 7 t)/{(p)t&a)- For the first term on the right 
hand side of Eq. (27) an additional minor assumption is made: 


_s_ ( avm\ „ _a_ ( B[F L /{p) t ) \ 

dxi \ f dxi j ~ Ox. y dxi J 


(28) 


This assumption is not necessary for the treatment of FMDF and is only adopted to establish 
consistency between the FMDF and the conventional LES. With these approximations, the 
modeled FMDF transport equation is 


m d[(u,) L F L ] 
dt dx{ 


_d_ 

dxi 


( 7 + it) 


dxi 


+ 


d^ c 


[D m (^» a — (^ct)l)Fl] — 


d[S a F L ) 

drpa 


(29) 


This equation may be integrated to obtain transport equations for the SGS moments. The 
equations for the first subgrid Favre moment, {4 >o,)l, and the generalized subgrid variance, 

= ( tf a )h ~ <0(«))l are 

djjpMM , d({p)j(uj) L (<l> a ) L ) _ _d_ 

dt dxi dx{ 


(7 + 7t) 


9 {^ a)l 

dii 


+ {p)t{S.)L 


(30) 


a(W ,<rl) + dUp)t(«i)LC,l) 


at 


dxi 


d_ 

dii 


(7 + It ) 


dal 


dxi 


+ 2(7 + It ) 


d{^(a))L d{<ft(a))l, 

dii dxi 


- 2 il m (p)tal + 2 (p)t (( <f>(c,)S{a))L - ( <f>(c))L{S( a ))L ) (31) 


where the subscripts in parenthesis are excluded from the summation convention. These 
equations are identical to those which can be derived by filtering Eq. (3) directly, and em- 
ploying consistent closures for the subgrid flux and the dissipation. In such direct moment 
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closure formulation, however, the terms involving (S a )z, remain unclosed. 

4 Monte Carlo Solution of the FMDF 

The Lagrangian Monte Carlo procedure (Pope, 1985) is employed for the solution of Eq. (29). 
In this procedure, each of the Monte Carlo elements (particles) obeys certain equations which 
govern their transport. These particles undergo motion in physical space by convection due 
to the filtered mean flow velocity and diffusion due to molecular and subgrid diffusivities. 
The compositional values of each particles are changed due to mixing and reaction. The 
spatial transport of the FMDF is represented by the general diffusion process governed by 
the stochastic differential equation (SDE) (Risken, 1989; Gardiner, 1990) 

dXi(t) = Di(X(t),t)dt + £(X(t), t)dWi(t) (32) 

where X is the Lagrangian position of a stochastic particle, £),- and E are the “drift” and 
“diffusion” coefficients, respectively, and Wi denotes the Wiener process (Karlin and Taylor, 
1981). The drift and diffusion coefficients are obtained by comparing the Fokker-Plank 
equation corresponding to Eq. (32) with the spatial derivative terms in the FMDF transport 
equation (Eq. (29)), 

E = ^2(7 + 7, )/(/>)/, + (33) 

The subgrid mixing and reaction terms axe implemented by altering the compositional 
makeup of the particles 

^ - (4.)i) + s„(<a + ) (34) 

where <f>+ = </> Q (X(t),t) denotes the scalar value of the particle with the Lagrangian position 
vector X{. The solutions of Eqs. (32) and (34) yield the same statistics as those obtained di- 
rectly from the solution of FMDF transport equation according to the principle of equivalent 
systems (Pope, 1985; Pope, 1994). 
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Numerical Solution Procedure 


A new computational algorithm is developed for the solution of the FMDF. While the al- 
gorithm is similar to that used in PDF methods (Pope, 1985), it is not exactly the same. 
Therefore, a detailed description is provided. 

The complete numerical solution of the equations governing the resolved field is based on 
a hybrid scheme in which the hydrodynamic Favre filtered equations (Eqs. (10)-(11)) are 
integrated by a finite difference method and the filtered scalar field is simulated by the Monte 
Carlo solution of the FMDF transport equation. The LES of the hydrodynamic variables, 
which also determines the subgrid viscosity and scalar diffusion coefficients, is conducted with 
the “compact parameter” scheme of Carpenter (1990). This scheme is based on a hyperbolic 
solver which considers a fully compressible flow. Here, the simulations are conducted at a 
low Mach number to minimize compressibility effects. All the finite difference operations 
are conducted on a fixed and uniform grid. Thus, the filtered values of the hydrodynamic 
variables are determined on these grid points. The transfer of information from these points 
to the location of the Monte Carlo particles (described below) is conducted via interpolation. 
Both fourth-order and second-order (bilinear) interpolations schemes were considered, but 
no significant differences in filtered quantities were observed. The results presented below 
utilize fourth- and second-order interpolation for two-dimensional (2D) and 3D simulations, 
respectively. 

The FMDF is represented by an ensemble of Monte Carlo particles, each with a set of scalars 
<f>^(t) = <fi a (X( n )(t),t) and Lagrangian position vector X^. A splitting operation is em- 
ployed in which transport in the physical and compositional domains are treated separately. 
The simplest means of simulating Eq. (32) is via the Euler-Maruyamma approximation (Kloe- 
den and Platen, 1995): Xl n \t k+l ) = X| n) (<*) -1- D\ n) (t k )At + El n \t k )(At) l ' 2 (\ n) {t k ), where 
At = t k+ i —t k is the computational time increment between two consecutive discretized time 
levels, Dj”^(f) = Di(X.^(t),t), E^(t) = E(X < ' n ^(t),t) and is a random variable with the 
standard Gaussian PDF. The coefficients D, and E require the input of the filtered mean 
velocity and the diffusivity (molecular and subgrid). These are provided by finite difference 
solution of Eqs. (10)-(11). 

The compositional values are subject to change due to SGS mixing and chemical reaction. 
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Equation (34) may be integrated numerically to simulate these effects simultaneously. Al- 
ternately, this equation is treated in a split manner. This provides an analytical expression 
for the subgrid mixing. Subsequently, the influence of chemical reaction is determined by 
evaluating the fine grained reaction rates and modifying the composition. 

The mixing model requires the Favre filtered scalar values. These and other higher moments 
of the FMDF at a given point are estimated by consideration of particles within a volume 
centered at the point of interest. Effectively, this finite volume constitutes an “ensemble 
domain” characterized by the length scale A E (not to be confused with A G ) in which the 
FMDF is discretely represented. A box of size A e is used to construct the statistics at the 
finite difference nodes. These are then interpolated to the particle positions. Since the SGS 
mixing model only requires the input of the filtered scalar values, and not their derivative, 
this volume averaging is sufficient. From a numerical standpoint, specification of the size of 
the ensemble domain is an important issue. Ideally, it is desired to obtain the statistics from 
the Monte Carlo solution when the size of sample domain is infinitely small (A E — ► 0) and 
the number of particles within this domain is infinitely large. With a finite number of parti- 
cles, if A# is small there may not be enough particles to construct reliable statistics. A larger 
ensemble domain decreases the statistical error, but increases the spatial error which mani- 
fests itself in artificially diffused statistical results. This compromise between the statistical 
accuracy and dispersive accuracy as pertaining to Lagrangian Monte Carlo schemes implies 
that the optimum magnitude of As cannot, in general, be specified a priori (Pope, 1985; 
Colucci et al., 1998). This does not diminish the capability of the scheme, but exemplifies 
the importance of the parameters which govern the statistics. 

In an attempt to reduce the computational overhead, a procedure involving the use of non- 
uniform weights is also considered. This procedure allows a smaller number of particles to 
be imposed in regions where a low degree of variability is expected. Conversely, in regions 
of highly varying character, a larger number of particles is allowed. This is akin to grid 
compression in finite difference (or finite volume) schemes. Operationally, the particles 
evolve with a discrete FMDF, 

N 

F N {i>- x, t) = Am £ - <£ (n) )6(x - x< n >) (35) 

n=l 

where u;f") is the weight of the n ttl particle and Am is the mass of a particle with unit weight. 
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The FMDF is the expectation of the discrete FMDF 


F L (V>;x,t) 


N 


Am ~ <£ (n) )6(x - x (n) )) 

n=l 

- <£ W )«(x - x<">)) 


( 36 ) 


for any n (1 < n < N). The brackets without the subscript L represent ensemble averaging. 
With integration of this expression over the composition domain within an infinitesimal 
volume, it is possible to demonstrate 


(p)t « 


Am 

AV 


E 

ti£Ae 


(37) 


where AV is the volume of the ensemble domain. The Favre filtered value of a transport 
quantity Q(<f>) is constructed from the weighted average 

£.«»< o'-W*’) 

1 ~ »<"> 

The approximations in Eqs. (37) - (38) are exact in the limit Ae — » 0 and the number 
of particles within the ensemble domain becomes infinite (Pope, 1985). Equation (37) im- 
plies that the filtered fluid density is directly proportional to the sum of the weights in the 
ensemble domain. With uniform weights, (p)t « and (Q)l ~ ^£Q(<£^) (Pope, 

1985) where Ne is the number of particles in the ensemble domain. Hence, with uniform 
weights, the particle number density decreases significantly in regions of high temperature. 
The implementation of variable weights allows the increase of the particle number density 
without a need to increase the number density outside of the reaction zone. 

To evaluate the chemical source terms, the fine grained values of the temperature (T^) for 
all particles are calculated from the composition variable <f> ^ = [V^' . . . , Y^J*\ 
and the fine grained values of density ( p (”)) are determined from evaluation of the equation 
of state at the reference pressure po. The filtered pressure is obtained by the filtered equation 
of state, (p)e = (p)e(TZT)i. In this equation ( p)( is obtained from the finite difference solver 
and the correlation (7 ZT)l is obtained by ensemble averaging in the Monte Carlo solver. In 
this way, the coupling between the hydrodynamic and the scalar fields is taken into account 
and allows the investigation of the effects of variable density. The results obtained by this 
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scheme axe identified by the label FMDF-1. 


The pressure ((p)/) field as determined by the above procedure exhibits some spatial oscil- 
lations caused by statistical error. Since the spatial derivatives of (p)t are required in the 
hydrodynamic solver, these oscillations can cause numerical difficulties. This is particularly 
exacerbated by the nature of the compressible hydrodynamic code which allows propagation 
of these oscillations throughout the computational domain. Our results shown below indicate 
that while the extent of noise in the pressure field is noticeable, it is not significant in the 
compositional variables. The amplitudes of the oscillations can be decreased by smoothing 
of the (TZT)l field. An alternate procedure is also followed in which the correlation (TZT)l 
is evaluated by the finite difference solution of its transport equation. With the assumption 
of constant TZ, only the solution of the Favre filtered temperature equation is required. The 
reaction source term in this equation is evaluated from the Monte Carlo solution. The re- 
sults obtained by this scheme are identified by the label FMDF-2. While the finite difference 
solution of the filtered temperature is used to calculate the filtered pressure in FMDF-2, 
the filtered temperature can also be evaluated directly from the Monte Carlo particles. The 
results below indicate that the filtered temperature fields obtained by the two methods axe 
nearly identical. 

In addition, another LES is also considered in which the modeled transport equations for the 
filtered scalar and the generalized SGS scalar variance are simulated with the finite difference 
scheme. The hydrodynamic solver and the models for the subgrid stress and mass flux axe 
identical to those in FMDF, but the effects of SGS fluctuations in the filtered reaction rate 
are ignored. That is, Eqs. (30)-(31) are solved via the finite difference scheme with the 
assumption ( S a (<f>))L — S a ((<f>)L)- The results based on this scheme are referred to as LES- 
FD. A variant of this model, in which the filtered reaction rate is modeled by ( S a (<f>))L — 
(5 0 ((</>)jr ))l was also considered. However this closure did not show any improvements 
over LES-FD; thus is not discussed. For non-reacting flows, the LES-FD results are used 
to demonstrate the consistency of the FMDF results. For reacting flows, the difference 
between FMDF and LES-FD demonstrates the effects of the SGS fluctuations. However, 
this comparison does not imply that these two methods are the only means of performing 
LES of reacting flows; several other schemes are currently available as indicated in Section 
1 . 


13 



Table 1: Attributes of the Computational Methods. 


Method 

Mean Field 
Equations 

Particle Properties 

Particle Fields 
Used in the Mean 
Field Equations 

Duplicate Fields 

LES-FD 

{p)t, («,)l, 

(<t>h 

— 

— 

— 

FMDF-1 

( p)t , («.)l 

<f>, P(<f> )> 

" 

(P)< 

FMDF-2 

(p)t, ( Ui ) L , 

(KT) l 

4>, 


< p )<, (ftru 


It is noted that the FMDF-1 simulation procedure is similar to that typically used in PDF 
methods (Pope, 1985; Tolpadi et al., 1995; Tolpadi et a/., 1996). The procedure described 
in FMDF-2 is proposed here for the first time. It is shown below that the pressure field 
as determined by this method exhibits almost no spatial oscillations, thus no smoothing is 
required. This scheme is starting to replace the equivalent of FMDF-1 in PDF methods 
(Pope, 1997). The attributes of the LES-FD, FMDF-1 and FMDF-2 schemes are outlined in 
Table 1. In this table, Snr denotes the source term in the equation governing the transport 
of nr. 


5 Results 

5.1 Flows Simulated 

The simulations of the following flow configurations are considered: 

1. A two-dimensional (2D) temporally developing mixing layer. 

2. A 3D temporally developing mixing layer. 

3. A 2D spatially developing planar jet. 

4. A 2D spatially developing mixing layer. 

The objectives of the numerical simulations are to: (i) demonstrate the consistency of the 
Monte Carlo solution procedure, (ii) demonstrate the capabilities of the FMDF, (iii) appraise 
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its overall performance, and (iv) highlight its deficiencies. The flow configurations (1),(3) 
and (4) are suitable for objectives (i) and (ii) in which 2D simulations are sufficient. How- 
ever, objectives (iii) requires 3D simulations. All flow configurations are used for objective 
(iv). The 2D simulations are conducted to allow extensive computations for assessing the 
consistency and accuracy of the FMDF and the convergence of the Monte Carlo results. 
Both non-reacting and reacting flows are simulated, and FMDF and LES-FD are applied 
to the cases itemized in (l)-(3). Some of these cases are also treated by DNS, the results 
of which are used to assess the performance of the FMDF. Further appraisal is made by 
comparison with laboratory data for the flow under item (4). 

The temporal mixing layer consists of two co-flowing streams traveling in opposite directions 
with the same speed (Riley et al., 1986; Jou and Riley, 1989; Givi, 1989). The reactants 
A and B are introduced into the top and the bottom streams, respectively. The length in 
the streamwise direction is large enough to allow for the roll-up of two large vortices and 
one (subsequent) pairing of these vortices. In 3D simulations, the length of the domain 
in spanwise direction is 60% of that in the streamwise direction. The layer is forced via 
both 2D and 3D forcing functions (Moser and Rogers, 1991; Miller et al . , 1994; Givi, 1994). 
The initial values of the reactants A and B at each spanwise location in 3D simulations are 
identical to those in 2D. In the figures presented below, x, y, z correspond to the streamwise, 
cross-stream and spanwise directions (in 3D), respectively in all the simulations. 

In the planar jet, the reactant A is issued from a jet of width D into a co-flowing stream 
with a lower velocity carrying reactant B (Givi and Riley, 1992; Steinberger et al., 1993). 
The size of the domain in the jet flow is 0 < x < 14Z), —3.5 D < y < 3.5 D. The ratio 
of the co-flowing stream velocity to that of the jet at the inlet is kept fixed at 0.5. A 
double-hyperbolic tangent profile is utilized to assign the velocity distribution at the inlet 
plane. The formation of the large scale coherent structures are expedited by imposing low 
amplitude perturbations at the inlet. The frequency of these perturbations correspond to 
the most unstable mode and subharmonics of this mode as determined by the linear stability 
analysis of spatially evolving disturbances (Michalke, 1965; Colucci, 1994). The characteristic 
boundary condition procedure developed by Poinsot and Lele (1992) is used at the inlet. 
This procedure facilitates evaluation of incoming waves which are necessary to satisfy the 
continuity equation. Zero derivative boundary conditions are used at the free-streams and 
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the pressure boundary condition of Rudy and Strikwerda (1980) is used at the outflow. 

The flow configuration in (4) is the one considered in the laboratory experiments of Mungal 
and Dimotakis (1984). In these experiments, a heat releasing reacting planar mixing layer 
composed of diatomic hydrogen in one stream and diatomic fluorine in the other stream is 
considered. Both reactants are diluted in nitrogen with the level of dilution determining the 
extent of heat release. While the laboratory flow, like all turbulent flows, is inherently 3D, 
it is dominated by large scale 2D structures (Brown and Roshko, 1974; Givi and Riley, 1992; 
Givi, 1994). We demonstrate that 2D simulations are sufficient to capture the hydrodynam- 
ics features of this flow reasonably well. The computational domain considers the region 
54.84 cm x 27.42 cm, which covers the whole region considered experimentally including 
x = 45.7 cm where measured data are reported. In order to mimic a “naturally” developing 
shear layer, a modified variant of the forcing procedure suggested by Sandham and Reynolds 
(1989) is utilized. The cross-stream velocity component at the inlet is forced at the most 
unstable mode as well as four harmonics (both sub- and super-) of this mode. A spatial lin- 
ear stability analysis was performed to determine the most unstable mode of the hyperbolic 
velocity profile imposed at the inlet. Sandham and Reynolds (1989) suggest the use of a 
random phase shift to “jitter” the layer and to prevent a periodic behavior. A similar ran- 
dom phase shift procedure is imposed here; a discrete approximation of the Wiener process 
is applied for the phase shift at each time increment. 

The flow variables are normalized with respect to selected reference quantities, denoted by 
the subscript r. In the temporal mixing layer, the reference quantities are the free-stream 
values and the reference length L r is defined such that ^ = 2.83, where £„ 0 is the initial 
vorticity thickness ( 6 V = r-° i — , where (tti)x, is the Reynolds averaged value of the 
Favre filtered streamwise velocity and A U is the velocity difference across the layer). In 
the spatial flows, normalization is performed with respect to the values in the high speed 
stream. In the planar jet L r — D. In the hydrogen-fluorine mixing layer, L r is equal to 
the distance from the virtual origin to the downstream measuring station in the experiment. 
These quantities are used to define the Reynolds number Re = SzMiLr p or the temporal 
mixing layer, the Reynolds number in terms of the total velocity difference across the layer 
(A U = 2 U T ) is Reg v0 = 5. 66 Re. The non-dimensional time is given by t* = ^. 
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5.2 Reaction Mechanisms 


For the flow configurations (l)-(3), the reaction scheme is of the type A + B — » V with 
an Arrhenius reactant conversion, = Sp = —pkfABexp(—E a /RT), where kj is the 
pre-exponential factor, E a is the activation energy, and A , B denote the mass fractions 
of species A, B, respectively. The species A, B, V are assumed thermodynamically iden- 
tical and the fluid is assumed to be calorically perfect. The normalized reaction rate is 
Sy — — p* DaAB exp(— Ze/T*) in which Ze = E a /RT r and Da = jpjfc denote the Zel- 
dovich number and the Damkohler number, respectively. T r denotes the reference ambient 
temperature. The degree of exothermicity is parameterized by the non-dimensional heat 
release parameter Ce = where A is the heat of reaction. Both constant rate and 

temperature dependent reactions are considered. 

The reaction mechanism associated with the mixing layer experiment is more complex. The 
hydrogen-fluorine reaction can be represented by the reaction (Mungal and Dimotakis, 1984) 

H 2 + F 2 -+ 2HF , A Q = 65 kcal/mol, (39) 

where A Q is the heat of reaction. This reaction belongs to the more general family of 
hydrogen-halogen reactions (Spalding and Stephenson, 1971; Chelliah, 1989). The heat 
released in a mixture containing 1% mole fraction of F 2 and 1% mole fraction of H 2 diluted 
in nitrogen results in an adiabatic temperature of 93 K above the ambient (Mungal and 
Dimotakis, 1984). The global representation in Eq. (39) is composed of a pair of second- 
order chain reactions (Mungal and Dimotakis, 1984) 

H 2 + F H F + H , AQ = 32 kcal/mol, i fc, = 2.6 x 10 12 T° 5 exp , (40) 

H + F 2 H F + F, AQ = 98 kcal/mol, k 2 = 3 x 10 9 T 15 expf ^ ■ j , (41) 

where the reaction rate constants k\ and k 2 are given in units of cm 3 / (mol s), T in K , and the 
universal gas constant R° in cal/(mol K). At low concentrations of the H atom, the reverse 
of the first of thesfe two reactions is slow. Additionally, the rate data suggest that the reverse 
of the second reaction is also negligible as compared to the forward reaction (Chelliah, 1989). 
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The explosion limits for the hydrogen-fluorine reaction indicate that a mixture of these two 
gases at typical ambient conditions is stable (Chen et al., 1975; Gmelin, 1980). Therefore, 
in order to initiate reaction, a source of F atoms must be provided (Mungal and Dimotakis, 
1984). Experimentally, this is accomplished by uniformly mixing a small amount of nitric 
oxide with the hydrogen-nitrogen mixture. The nitric oxide reacts with the fluorine to 
produce free fluorine atoms 

NO + F 2 NOF + F, A Q = 18 kcal/mol, k 3 = 4.2 x 10 11 exp (~ ^) * ( 42 ) 

The reverse of this reaction may be neglected (Rapp and Johnston, 1960). An additional 
reaction serves to limit the nitric oxide concentration (Baulch et al., 1981; Cool et al., 1970) 

F + NO + M -^4 NOF + M, AQ = 57 kcal/mol, k 4 « 3 x 10 16 cm 6 /(mol 2 s). (43) 

While it is necessary to add nitric oxide to initiate reaction, the addition of excessive amounts 
would deplete the availability of the free F atoms. Mungal and Dimotakis (1984) indicate 
that keeping the product of nitric oxide and diatomic fluorine molar concentrations at 0.03% 
results in a rapid combustion. It was also noted that an increase of 50% in the nitric oxide 
concentration results in no appreciable changes in the temperature. This suggests that the 
hydrogen-fluorine reaction can be approximated by the limit of infinite rate chemistry. In 
the simulations, therefore, both finite and infinite rate models are considered. Due to the 
very fast rate of the reaction, the compositional change due to reaction is implemented 
in 10 incremental time steps for every hydrodynamic time step. These simulations with 
stiff reaction rates are obviously computationally intensive. The implementation of the 
infinite rate chemistry model (Williams, 1985) is significantly less expensive. With this 
approximation, it may be possible to employ the assumed FDF approach (Madnia and Givi, 
1993). However, in order to demonstrate the operationally of the FMDF, here this procedure 
is employed for both finite and infinite rate models. 

5.3 Numerical Specifications 

The magnitude of the flow parameters considered in DNS are dictated by the resolution 
which can be afforded. The primary parameters are Re, Da, Ze, Ce, Sc, and Pr. In 
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all simulations Sc = Pr = 1. All finite difference simulations (in both DNS and LES) are 
conducted on equally-spaced grid points (Ax = Ay = A z (for 3D) = A). The highest 
resolution in DNS of the 2D temporal mixing layer consists of 433 x 577 grid points which 
allows reliable calculations at Re = 2,000, Ce = 5, Ze = 8, and Da = 11.92. The DNS of 
the 3D temporal shear layer is conducted with a resolution of 217 x 289 x 133 grid points 
with Re — 400, Da = 1 and Ce = Ze — 0. The DNS of the planar jet is performed on 
1201 x 601 grid points and allows accurate simulations with Re — 10,000, Ce = 2.5, Ze = 8 
and Da = 119.2. The FMDF and LES-FD axe conducted with lower grid resolutions. The 
LES of the temporal mixing layer is conducted on 37 x 49 and 55 x 73 grid points for 2D 
simulations while resolutions of 37 x 49 x 23 and 55 x 73 x 34 are utilized in 3D. The LES 
of the spatial jet and hydrogen-fluorine mixing layer are conducted on 201 x 101 grid points. 
A top-hat filter function (Aldama, 1990) of the form 
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is used with Ac = 2A and Nd denotes the number of dimensions. No attempt is made to 
investigate the sensitivity of the results to the filter function (Vreman et al., 1994) or the 
filter size (Erlebacher et al ., 1992). 


For FMDF simulations of the temporal mixing layer, the Monte Caxlo particles are initially 
distributed throughout the computational region. For the jet flow, the particles are supplied 
in the inlet region — 1.75D < y < 1.75 D. As the particles convect downstream, this zone 
distorts as it conforms to the flow as determined by the hydrodynamic field. In regions pop- 
ulated with particles, 5ZneA E remains proportional to the instantaneous filtered density 
(within statistical error). In regions without particles, a delta function FDF corresponding to 
the free-stream composition is enforced. The simulation results are monitored to ensure the 
particles fully encompass and extend well beyond regions of non-zero vorticity and reaction. 
In the temporal mixing layer, due to flow periodicity in the streamwise and spanwise direc- 
tions, if the particle leaves the domain at one of the boundaries new particles are introduced 
at the other boundary with the same compositional values. In the spatially evolving jet and 
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the planar mixing layer, new particles are introduced at the inlet at a rate corresponding to 
the desired (imposed) local particle number density and fluid velocity. In some of the planar 
jet simulations and all of the hydrogen-fluorine mixing layer simulations, variable particle 
weights are employed. With prescription of the filtered fluid density, the particle weight is 
adjusted to yield the proper mass flux across the boundary. All other simulations utilize 
uniform weights. The sensitivity of the statistical results to the number of particles per grid 
cell ( NPG ) and the size of the ensemble domain is studied in the temporal mixing layer. 
The following sizes are considered: A e = 2A, A, A/2. 

In the FMDF simulation of the experimental mixing layer configuration, initially N PG = 5 
in the free-streams and gradually increases in the cross steam direction peaking to NPG = 25 
at the splitter plate (y = 0). This yields 20 to 100 sample points per ensemble for As = 2A. 
The particles are supplied in the region — 0.12T r < y < 0.12T r where L r = 45.7 cm. The 
composition of incoming particles is set according to composition of the fluid at the point of 
entry. The magnitudes of the Reynolds, Peclet, Damkohler and Zeldovich numbers and the 
velocity ratio across the layer in the simulations are the same as those in the experiment, 
but the maximum value of the Mach number in the simulations is 0.31 which is higher than 
that in the experiment. This was necessary in the compressible flow solver employed for 
the simulations. With the values of the physical parameters in this experiment, it is not 
possible to employ DNS and LES-FD for this flow, thus only FMDF results are compared 
with experimented data. For that FMDF-1 is used in which smoothing of (7 ZT)l is done 
with a box filter consisting of 3 x 3 grid points with equal weights. 

The simulated results are analyzed both instantaneously and statistically. In the former, the 
instantaneous contours (snap-shots) and the scatter plots of the scalar values are considered. 
In the latter, the “Reynolds-averaged” statistics are constructed form the instantaneous data. 
In the temporal mixing layer, the statistics are constructed by the ensemble from all the grid 
points in the homogeneous direction(s) x (and 2 in 3D). In the spatially developing mixing 
layer and the jet flow, averaging is conducted via time sampling. All Reynolds averaged 
results are denoted by an overbar. In the presentations below, the asterisk (denoting the 
normalized quantities) is dropped. 
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5.4 Model Parameters 


In the implementation of the MKEV, the magnitude of the reference velocity U, is set to zero 
in the cross-stream and spanwise directions, and to the average of the high and low speed 
streams in the streamwise direction. Additionally, the ratio of the filter size at the secondary 
level to that at the grid level is A g>I&g = 3. In all simulations (7/ = 0.006. The magnitude 
of Cr is 0.020 and 0.013 for 2D and 3D, respectively. The subgrid mass flux is modeled via 
Eq. (14). In all cases except LES of the hydrogen-fluorine mixing layer, Pr t = Set = 0.7. No 
attempt is made to determine the magnitudes of these model constants in a dynamic manner 
(Germano, 1992). The subgrid mixing model requires the input of the constant Cq which 
also determines the SGS variances. The value Cq = 4 is used in most simulations. In the 
hydrogen-fluorine configuration Set — Pf't — 0-4 and Cq = 6. Some constant density test 
simulations are also conducted in which Cq = 3 as previously used by Colucci et al. (1998). 
The non-universality (flow dependence) of the hydrodynamic model constants (C/, Cr, Pr t , 
Set) has been well recognized and was expected here. The additional constant introduced 
by FMDF is Cq, although this. constant also appears if the SGS variance is considered in the 
conventional LES-FD. This non- universality, in general, diminishes the predictive capability 
of LES; however the range of the values as considered here is not very broad. 

5.5 Consistency of FMDF 

The objective in the results presented in this subsection is to demonstrate the consistency 
of the FMDF formulation. For this purpose, the LES results via FMDF and LES-FD are 
compared against each other in 2D and 3D temporal mixing layers. Since the accuracy of 
the finite difference scheme is well-established, this comparative analysis provides a means of 
assessing the performance of the Monte Carlo solution of the FMDF. For most of the results 
in this section, NPG = 50 in 2D and NPG = 20 in 3D at locations where ( p)t = 1. In 
2D, A e = A and in 3D, A# = 2 A. Several additional simulations are also performed with 
varying values of NPG and A e to asses their effects. 

Simulations of 2D non-reacting temporally developing mixing layers are conducted in which 
the flow is initiated with non-uniform density and temperature distributions. The initial 
filtered density is distributed as a “spike” at the middle of the layer. With uniform weights 
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assigned to the Monte Carlo particles, the particle number density must remain proportional 
to the fluid density. This is observed in Fig. 1, where it is demonstrated that the filtered 
density evaluated from the Monte Carlo particles matches very well with that of the finite 
difference calculated values at the Eulerian grid points. The values generated by the finite 
difference solution are identified by FD and the results generated by ensemble averaging of 
the Monte Carlo particles are identified by MC. Figure 1 shows that at the final time of 
the simulation (when the flow has experienced the pairing of two neighboring vortices) the 
Reynolds averaged filtered density calculated by the finite difference and the Monte Carlo 
procedures are very close. The particle number density exhibits an appreciable degree of 
oscillations due to statistical errors associated with a finite sample of particles. 

Figure 2 shows the temporal evolution of the vorticity thickness. When the flow starts with 
uniform density, the effect of thermodynamic quantities on the hydrodynamics is negligible. 
Thus the S v profiles as obtained by FMDF-1 and FMDF-2 are nearly identical. With an 
initial density spike, the growth of the layer is damped as expected (McMurtry et al., 1989; 
Jackson, 1992; Colucci, 1994), but the results obtained by FMDF-1 are very close to those 
by FMDF-2. The slight differences are due to the numerical solution procedures. The results 
obtained by both procedures are close to those obtained by DNS. 

In Fig. 3, the contour plots of the resolved vorticity and temperature at the fined time 
(t = 44) as obtained by FMDF-1 and FMDF-2 are shown. This figure provides a visual 
demonstration of the consistency of the FMDF as the results via the two FMDF procedures 
are similar. The difference, as expected, is exhibited by the oscillations in FMDF-1. The 
effect of the baroclinic torque in generating vorticity near the braids is captured by both 
simulations. To exhibit the extent of the noise more clearly, the Reynolds averaged values of 
the resolved pressure and the mass fraction of a conserved scalar are shown in Fig. 4. The 
most significant difference is evident in the filtered pressure field which exhibits appreciable 
oscillations in FMDF-1. These oscillations are reduced by application of a local least square 
filter to smooth the Monte Carlo ( T)l field. This operation does not modify the other 
statistical quantities. Several other filter functions are also considered and their influence 
is summarized in Fig. 5 where the percentages of the differences between the values of (p)/ 
via FMDF-2 and FMDF-1 with smoothing are shown. In all cases, the difference is small 
(less than 2%); the most significant difference is expectedly observed when no smoothing 
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operation is applied. Figure 5 also shows that the difference is significantly decreased as the 
number of Monte Carlo particles is increased. 

To demonstrate the consistency between the FMDF and LES-FD, a comparison is made 
between the moments of the mass fraction of A in the non-reacting temporal mixing layer 
with an initial density spike as obtained by the two procedures. Figure 6 shows the instan- 
taneous contour plots of the Favre filtered mass fraction of species A and Fig. 7 shows the 
Reynolds averaged values of the moments of this mass fraction. In these simulations, the 
filtered temperature is calculated via FMDF-1 without smoothing. The similarity of FMDF 
and LES-FD results is evident in both figures. The agreement in the first moment (Figs. 
7(a,c)) is quite good even for large values of A# and small values of NPG. The difference is 
more apparent in the subgrid variance values (Figs. 7(b,d)). However, the difference becomes 
smaller as A e decreases. 

In reactive flows, the consistency established above no longer exists since the reaction term 
appears in a closed form in the FMDF formulation but not in the moment equations of LES- 
FD. This inconsistency, which motivates the use of FMDF, is illustrated in Fig. 8 where the 
temporal evolution of the integrated “total product” (fip(t) = f (P)l(j/, t)dy^ in a constant 
density reacting temporal mixing layer with Da = 2 and Ce — Ze = 0 is shown. In these 
simulations, the LES resolution is 37 x 49 and Re = 500. The LES results are also compared 
with those obtained via DNS with 433 x 577 grid points. It is shown that the FMDF 
results are very close to those via DNS, but LES-FD significantly overpredicts the amount 
of products formed. Also shown in Fig. 8 are the results via the constant density filtered 
density function (FDF) formulation (Colucci et al., 1998) which is suitable for this flow. The 
close agreement of FMDF, FDF and DNS results indicate both the consistency of the Monte 
Carlo solution and the relative superiority of FMDF over LES-FD. 

To generalize the results above, LES of a 3D temporally developing mixing layer is con- 
ducted. In these simulations, a non-reacting flow with a density spike similar to that in 
2D is considered. The statistical results in simulations with 3D forcing exhibit significant 
variations along the spanwise direction. The filtered pressure obtained from FMDF-1 ex- 
hibits similar trends to those obtained from FMDF-2 but does portray statistical noise. As 
is the case for 2D simulations, the filtered mass fraction and temperature calculated by the 
Monte Carlo solver are close to those obtained by the finite difference simulations. This is 
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illustrated Fig. 9 in which scatter plots of { T)l and (A)l values generated by FMDF-2 are 
shown. The correlation coefficient between the Monte Carlo (MC) and the finite difference 
(FD) values is 0.999 for both sets of results shown. 

5.6 Validation via DNS Data 

The objectives in this subsection are to assess the overall performance of FMDF, to appraise 
the validity of the submodels employed in the FMDF transport equation, and to demonstrate 
the capabilities of FMDF for LES of exothermic chemically reacting flows. To meet these 
objectives, the FMDF results are compared against DNS results of the same flow configura- 
tions with the same magnitudes of the physical parameters (Re, Da, etc.). For a meaningful 
comparison, the DNS data are filtered and down-sampled onto coarse grid points correspond- 
ing to those employed in FMDF. At this point it is emphasized that FMDF is not claimed 
to be an alternative to DNS; the comparisons made here are primarily for assessment of 
the FMDF. For further comparative assessments, the FMDF results are also compared with 
those via LES-FD. Both 2D and 3D simulations are considered. Unless otherwise specified, 
all Monte Carlo simulations presented in this section are based on the FMDF-2 formulation. 

To quantify the performance of FMDF in LES of the exothermic reacting 2D temporal mixing 
layer, in Fig. 10 the cross-stream variation of the Reynolds averaged filtered temperature 
values at t = 44 are shown. In this simulation, Da = 11.92, Ze = 8 and Ce = 5. The FMDF 
results are calculated with both A e = A and A e =■ 2 A. Initially, the particle number density 
is set to N PG = 40 with initial uniform fluid density. The size of the ensemble domain for 
the evaluation of the Favre filtered statistics does not have a significant influence on the first 
filtered moment. The deviation of LES-FD results from those via FMDF and/or DNS is 
evident. This behavior is observed at all times for all the cases considered. It is expected 
that the difference between DNS and LES-FD results would be even more as the magnitude 
of the Damkohler number and/or Reynolds number increases (Colucci et al., 1998). Figure 
10 shows that for this flow with a rather significant variation of temperature, the averaged 
filtered temperature is predicted well by FMDF. Comparatively, LES-FD overpredicts the 
filtered temperature values. While the finite difference solution of the filtered temperature 
is used to calculate the filtered pressure in the FMDF-2, the filtered temperature can also 
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be evaluated directly from the ensemble of the Monte Carlo elements. Figure 10 indicates 
that the evaluation of the filtered temperature in this way (denoted by MC ensemble) is 
consistent with that obtained by FMDF-2. 

The results of the spatially developing jet flow are shown in Figs. 11-17 in which severed 
issues pertaining to the Monte Carlo simulation are addressed. Figures 11 and 12 show the 
instantaneous contours of the filtered pressure and the filtered temperature values, respec- 
tively. Parts (a), (b), and (c) of these figures correspond to results with FMDF-1 without 
smoothing of the temperature field, FMDF-1 with smoothing, and FMDF-2, respectively. 
While the temperature fields as obtained by all three procedures are similar, the differences 
between the pressure fields axe noticeable. The behavior portrayed in Fig. 11(c) is physical, 
whereas the oscillations observed in Figs. ll(a,b) could cause numerical problems. While 
these oscillations did not cause problems here, Fig. 11 shows that FMDF-2 is more robust 
and is recommended for both LES and PDF simulations. In Figs. 13 and 14 the influence 
of the particle weights in the Monte Carlo simulation is exhibited. Figure 13 shows that the 
instantaneous particle number density and the filtered fluid density calculated by FMDF are 
highly correlated in these simulations in which uniform particle weights are employed. It 
is noted that the particle number density is lowest in the high temperature reaction zones. 
Figure 14 shows the results via variable weights. It is observed in Fig. 14(a) that there 
is a higher concentration of particles in the reaction zones in comparison to the case with 
uniform weights. The particle mass density shown in Fig. 14(b) is highly correlated with the 
filtered fluid density (Fig. 14(c)). A comparison between Figs. 13(b) and 14(c) indicates that 
despite the significant difference in the total number of particles and particle weighing pro- 
cedures, the filtered density fields are nearly identical in the two simulations. This similarity 
is also reflected in the streamwise variations of the total product [Sp{x ) = / (P)l(x, y)dy^j 
in Fig. 15. The results via both procedures are nearly identical and are superior to LES-FD 
in matching with DNS results. The computational time in the simulations with variable 
weights is about half of that in simulations with equal particle weights. 

As indicated previously, the essential difference between FMDF and LES-FD is due to the 
ability of FMDF in accounting for the SGS scalar fluctuations. To demonstrate this explicitly, 
in Fig. 16, the contour plots of the “SGS unmixedness” defined as ( p)t [(£(<£))!, — <S((<£)t)] 
are shown. It is observed that the FMDF results are in good agreement with DNS. The con- 
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tribution of the SGS unmixedness to the toted filtered reaction rate is expected to increase as 
the magnitudes of the Re , Da , Ce increase. Therefore, it is anticipated that the difference 
between DNS and LES-FD results would be even more with increased values of these param- 
eters. Scatter data of the instantaneous product mass fraction V vs. the mixture fraction 
Z are presented in Fig. 17. These data are gathered at the final time of the simulations 
including the results within the region x > 3.5J9. Both the DNS (Fig. 17(a)) and FMDF-2 
(Fig. 17(b)) exhibit significant scatters indicative of appreciable finite rate chemistry effects. 
The FMDF is able to capture the scatter reasonably well. It is important to note that while 
the fine-grained values associated with the particles may be interpreted as instantaneous 
realizations, conventional LES cannot offer such “de-filtered” information. 

The major conclusions drawn from the 2D results axe confirmed in 3D simulations. In Fig. 
18, the time-variation of the total product as predicted by FMDF of the constant density 
temporally developing reacting mixing layer is compared with DNS and LES-FD results. 
Consistent with the 2D results, the total product predicted by FMDF is closer to DNS in 
comparison to that of LES-FD. With increased resolution in LES, the difference between 
DNS and LES-FD is less, but the FMDF results are not significantly modified. 

5.7 Validation via Laboratory Data 

The experiments of Mungal and Dimotakis (1984) were conducted with several values of the 
equivalence ratio, <f> = C 02 /C 01 where Co refers to the free-stream molar concentration and the 
subscripts 1 and 2 denote the reactants in the high- and the low-speed streams, respectively. 
Equivalence ratios of 1, 2, 4 and 8 were considered. In addition, “flip” experiments were also 
conducted in which inverse values of the equivalence ratio (<f> = 1, |, ^ and |) were consid- 
ered. All of these cases are considered in the simulations by FMDF-1. The implementations 
of DNS and LES-FD axe not possible for this flow. 

Figure 19 displays the contour plots of the instantaneous and the Reynolds averaged tem- 
perature field for the case with 4> = 1. In this simulation, the finite rate reaction scheme is 
employed. The peak value of the instantaneous temperature field approaches, but is lower 
than, the adiabatic flame temperature. This is due to the filtering of the temperature field. 
The peak values of the time averaged temperature values are considerably lower than that 
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of the adiabatic flame temperature, an intuitive fact indicated by Mungal and Dimotakis 
(1984) and also by Wallace (1981). However, it is noted that a large number of individual 
particles (i.e. realizations) do indeed approach the adiabatic limit. 

The FMDF predictions are compared with experimental results both qualitatively and quan- 
titatively. Figure 20 shows the time history of the temperature at several cross stream loca- 
tions as obtained by FMDF. Each vertical increment represents temperature values ranging 
from the ambient to the maximum attained instantaneous temperature (T max ). These time 
traces are quaditatively similar to those measured experimentally (Mungal and Dimotakis, 
1984). One notable difference is observed near the middle region of the layer. In this region, 
there are instances when the simulations exhibit near ambient temperature vadues (cold 
fluid). While there is some evidence of this behavior in the experiments, it is more pro- 
nounced in the simulations. This is partly attributed to the 2D nature of the simulations as 
the small scale mixing present in 3D tend to provide a more effective mixing (Miller et a/., 
1994). For this reason it is expected that the minimum values of the time averaged temper- 
ature in the vicinity y = 0 to be slightly lower than those measured experimentadly. Another 
reason for this difference is due to the fact that the cold wire probes may include some 
thermal lag and conduction errors (Scadron and Warshawsky, 1952; Paranthoen et al., 1982; 
Mungal and Dimotakis, 1984) manifesting in an artificial “smoothing” effect in the measured 
temperature values. 

For a quantitative comparison, in Fig. 21 the cross stream variations of the Reynolds averaged 
temperature rise normalized by the adiabatic temperature rise (T 0 ) are shown. The quantity 
S\ denotes the distance between the points where the cross stream mean temperature rise is 
1% of the maximum mean temperature rise and j/o is the cross stream location where the 
time-averaged streamwise velocity is the average of the high- and low-speed velocities. No 
attempt is made to de-filter the LES results and ( T)l is directly compared to experimented 
data. The agreement between the FMDF and experimental data is good. Also shown in this 
figure are the results based on the FMDF with the infinite reaction rate model. As expected, 
the results are very close to those of the finite rate simulation, but the computational cost 
is significantly less. In this particular case, the time requirement for FMDF simulations 
with the infinite rate chemistry is approximately 16% of that for the finite rate chemistry 
simulations. Due to this lower cost, and the confidence in the infinite rate model, the 
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remaining simulations are conducted with this model. 

To demonstrate the flip effect, Fig. 22 shows the cross-stream variation of the normalized 
temperature for all equivalence ratios (the simulations with <f> = 1 are repeated). Two obser- 
vations are made consistent with the experimental results: (1) the peak value of the mean 
temperature in each of the experiments is different from that in the corresponding flip ex- 
periment, although the adiabatic flame temperature is the same, (2) the peak temperature 
value shifts toward the lean reactant stream. Since the only difference between each of the 
two cases is the interchange of the low and high speed reactants, the reason for this behavior 
is due to the different entrainment processes (Mungal and Dimotakis, 1984). Additionally, 
with the exception of the two cases with <f> = 1, the peak temperature is higher for equiva- 
lence ratios greater than one compared to the reciprocal equivalence ratios. Consistent with 
the experimental results, the peak normalized temperature reaches a maximum for an equiv- 
alence ratio in the range 1 < <j> < 2. These trends are more clearly portrayed in Fig. 23(a), 
which Mungal and Dimotakis (1984) refer to as “inferred” temperature profiles. These reflect 
the temperature if the high speed reactant was fixed at 1 % molar concentration while the 
low speed stream was varied from |% through 8% to obtain the desired equivalence ratios. 
This figure supports the conclusion of Mungal and Dimotakis (1984) that there exists an 
asymptotic limit to the amount of products formed as the high speed reactant is burned to 
completion. A similar behavior is exhibited in Fig. 23(b) in which the inferred temperature 
profiles are shown for the situation in which the low stream reactant is fixed at 1% and the 
high speed reactant is varied to obtain the same equivalence ratios. 

Further quantitative comparison between the FMDF and experimental results is made in 
Fig. 24 which shows the variation of the normalized product thickness with the equivalence 
ratio. The product thicknesses are defined as (Mungal and Dimotakis, 1984) 
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where C v is the molar heat capacity of the carrier gas. Figure 24(a) indicates the FMDF 
predicts the extent of product formation reasonably well over a wide range of equivalence 
ratios. At low values of <f > , the amount of products varies nearly linearly as the low speed 
reactant is consumed when excessive amounts of the high speed reactant are present. At 
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Table 2: The computational times for 2D planar jet simulations. 


Simulation 

Grid resolution 

Normalized CPU timef 

DNS 

1201 x 601 

242.5 

FMDF 

1 201 x 101 

7.62 

LES-FD 

201 x 101 

1 


f Unit correspond to 760 seconds on a Cray-C90. 


high values of the equivalence ratio, the product thickness approaches an asymptotic value 
as the reaction progress is inhibited by a lack of high speed reactant relative to the amount 
of reactant in the low speed stream. Figure 24(b) demonstrates a similar agreement between 
the experimental and the FMDF results. 


5.8 Computational Requirements 

To appraise the computational requirements of the FMDF, the computational times for some 
of the cases are monitored. Tables 2 and 3 list the normalized CPU times required for the 
simulations of the reacting 2D planar jet and the reacting 3D temporally developing mixing 
layer, respectively. These cases are selected since simulations via all three schemes (FMDF, 
LES-FD and DNS) are conducted. The computational times listed for FMDF are those asso- 
ciated with FMDF-2, although the increase in cost over FMDF-1 is insignificant. Obviously 
the overhead of the FMDF simulation is extensive as compared to LES-FD; nevertheless, the 
computational time for FMDF simulation is significantly less than that of DNS. Again it is 
emphasized that FMDF is not claimed to be an alternative to DNS; neither it is claimed that 
the FMDF is capable of reproducing all DNS results. However, the close proximity of values 
obtained FMDF and DNS, and the substantial lower computational costs of FMDF makes 
it as a viable tool for simulations of reacting flow systems for which DNS is not possible. 


6 Summary and Concluding Remarks 

The basic objective of this work is to develop a methodology for large eddy simulation (LES) 
of turbulent reacting flows, with inclusion of exothermicity and variable density effects. The 
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Table 3: The computational times for 3D temporal mixing layer simulations. 


Simulation 

Grid resolution 

Normalized CPU timef 

DNS 

217 x 289 x 133 

182.71 

FMDF 

55 x 73 x 34 

7.64 

LES-FD 

55 x 73 x 34 

1 


f Unit correspond to 655 seconds on a Cray-C90. 

methodology is termed the “filtered mass density function” (FMDF) and is based on the 
extension of the “filtered density function” (FDF) developed previously for LES of constant 
density, reacting, isothermal flows (Colucci et al ., 1998). The procedure for this extension is 
similar to that used in probability density function (PDF) methods in Reynolds averaging 
procedures (Pope, 1985). Here the FMDF is considered for treatment of scalar variables. A 
transport equation is developed for the FMDF in which the unclosed terms, similar to PDF 
methods, are due to SGS convection and mixing. The former is modeled via the gradient 
diffusion model as done in most LES of non-reacting flows (Galperin and Orszag, 1993); the 
latter is closed via the IEM model as typically used in PDF methods (Pope, 1985). 

The modeled FMDF transport equation is solved numerically via a Lagrangian Monte Carlo 
scheme in which the solutions of the equivalent stochastic differential equations (SDEs) are 
obtained. Two Monte Carlo procedures are considered. The first (FMDF-1) is similar to that 
typically used in PDF methods (Pope, 1985; Tolpadi et al., 1995; Tolpadi et al., 1996). The 
second (FMDF-2) is new. Both schemes preserve the Ito-Gikhman nature of the SDEs and 
provide a reliable solution for the FMDF. The second scheme is more robust in dealing with 
the statistical noise generated by the Monte Carlo scheme. The consistency of the FMDF, 
the convergence of its Monte Carlo solutions, the advantages and drawbacks of the FMDF 
as well as the performance of the closures employed in the FMDF transport equation tire 
assessed. This is done via extensive comparisons between the results obtained by the Monte 
Carlo procedure and the finite difference solution of the transport equations of the first two 
filtered moments of scalar quantities (LES-FD). In non- reacting flows, the consistency and 
convergence of the Monte Carlo solution is demonstrated by good agreements of the first two 
SGS scalar moments with those obtained by LES-FD. The performance of FMDF and its 
superiority over LES-FD are demonstrated by comparison with direct numerical simulations 
(DNS) results of two-dimensional (2D) and 3D temporally developing mixing layers, and 
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a 2D spatially developing jet. In all cases the FMDF results are shown to be in closer 
agreement with the DNS data than are the LES-FD results in which the influence of the 
SGS fluctuations on the reaction rate are ignored. 

The performance of the FMDF is further appraised by comparison against the experimental 
data of Mungal and Dimotakis (1984) of a spatially developing mixing layer involving the 
exothermic hydrogen-fluorine reaction. The FMDF is considered via both finite rate and 
infinitely fast chemistry. The treatment of the former with a stiff reaction source term 
is computationally expensive, but comparison of the results with those of the latter gives 
confidence in the less costly infinite rate procedure. The results produced by both methods 
compare favorably with experimental data find some qualitative features, such as the “flip 
effect”, are captured by the FMDF simulation. 

In addition to the those in the hydrodynamic closure, there axe three constants for the 
LES of scalar quantities: Set and Pr t for the SGS convective fluxes of the mass fraction 
and the temperature, respectively, and Co as appears in the SGS mixing model. Based 
on the present results and those of Colucci et al. (1998) for a variety of different flows 
(2D and 3D, constant and variable density, different chemistry schemes, etc.) it seems that 
Set = Prt ~ 0.4 — 0.7, Cn « 3 — 6. The predictive capability of the FMDF can be improved 
by future developments in PDF methods. 

While the FMDF method is computationally more expensive than conventional LES method, 
it is much more advantageous for treating reacting flows. The computational overhead is 
tolerable for simulations of complex reacting flows for which DNS is not feasible. 
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Figure Captions 


Figure 1. Cross-stream variation of the filtered density in the 2D temporal mixing layer 
obtained by FMDF-1 at t = 44. 

Figure 2. Vorticity thickness vs. time in the 2D temporal mixing layer. 

Figure 3. Contours of the filtered vorticity and temperature in the 2D temporal mixing layer 
obtained by FMDF-1 (right side) and FMDF-2 (left side) at t = 44. Top: vorticity field, 
bottom: temperature field. 

Figure 4. Cross-stream variation of the mean filtered quantities at t = 44 in the 2D temporal 
mixing layer, (a) Pressure, (b) mass fraction of a conserved scalar. 

Figure 5. Cross-stream variation of the percentage of the difference in pressure as obtained by 
FMDF-2 and FMDF-1 with different smoothing in the 2D temporal mixing layer at t = 44. 
Long-dashed line: no smoothing, NPG = 50 and A# = A; Dotted-dashed line: smoothed 
with a Gaussian filter, NPG = 50 and Ae = A; solid line: smoothed with a box filter, 
NPG = 50 and Ae = A; dashed line: smoothed with a local least square filter, NPG = 50 
and A e — A; Long-dashed thick line: smoothed with a Gaussian filter, NPG = 200 and 
A e = A; Dotted thick line: smoothed with a Gaussian filter, NPG = 50 and A e — 2A. 

Figure 6. Contours of the filtered values of the conserved scalar at t = 44 in the 2D temporal 
mixing layer as obtained by (a) LES-FD, (b) FMDF. 

Figure 7. Cross-stream variation of mean filtered scalar ((a) and (c)) and the generalized 
variance of the conserved scalar ((b) and (d)) in the 2D temporal mixing layer. 

Figure 8. Total product variation with time in the 2D temporal mixing layer. 

Figure 9. Scatter plots of the filtered quantities as obtained by the Monte Carlo (MC) 
solution vs. those via the finite difference (FD) solution in the 3D temporal mixing layer: 
(a) temperature, (b) the conserved mass fraction. 

Figure 10. Cross-stream variation of the normalized filtered temperature in the 2D temporal 
mixing layer at t = 44. 

Figure 11. Contours of the normalized filtered pressure in the planar jet: (a) FMDF-1 with 
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no smoothing of the filtered temperature, (b) FMDF-1 with smoothed filtered temperature 
with box filter, (c) FMDF-2. 

Figure 12. Contours of the normalized filtered temperature in the reactive planar jet: (a) 
FMDF-1 with no smoothing of the filtered temperature, (b) FMDF-1 with smoothed filtered 
temperature with box filter, (c) FMDF-2. 

Figure 13. Contours of (a) the particle number density, (b) the fluid filtered density in the 
reactive planar jet simulations with uniform weights. 

Figure 14. Contours of (a) the particle number density, (b) the particle mass density, (c) the 
fluid filtered density in the reactive planar jet simulations with variable weights. 

Figure 15. Streamwise variation of the total product in the reactive planar jet. 

Figure 16. Contours of the normalized instantaneous SGS unmixedness in the reactive planar 
jet, (a) DNS, (b) FMDF. 

Figure 17. Scatter plot of the product mass fraction vs. the mixture fraction in the reactive 
planar jet, (a) DNS, (b) FMDF. 

Figure 18. Total product vs. time in the 3D temporal mixing layer, (a) lower LES resolution 
(37 x 49 x 23), (b) higher LES resolution (55 x 73 x 34). 

Figure 19. Contour plots of (a) instantaneous Favre filtered temperature, (b) time averaged 
Favre filtered temperature for (j) = 1 in the hydrogen-fluorine mixing layer. The values are 
normalized by T r . 

Figure 20. Time history of the instantaneous Favre filtered temperature in the hydrogen- 
fluorine mixing layer at several cross stream locations. 

Figure 21. Cross stream variation of the normalized mean temperature for <f> = 1 in the 
hydrogen-fluorine mixing layer. 

Figure 22. Cross stream variation of the normalized mean temperature for all equivalence 
ratios in the hydrogen-fluorine mixing layer. 

Figure 23. Cross stream variation of the “inferred” mean temperature profiles for (a) 1% 
high speed mole fraction, (b) 1% low speed mole fraction for all equivalence ratios in the 
hydrogen-fluorine mixing layer. 
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Figure 24. Normalized product thickness variation with equivalence ratio in the hydrogen- 
fluorine mixing layer: (a) <5 pl vs. the equivalence ratio, (b) S p2 vs. the inverse equivalence 
ratio. 
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1. Introduction 

Large eddy simulation (LES) of turbulent reacting flows has been the 
subject of widespread investigation (McMurtry et ai , 1992; Galperin and 
Orszag, 1993; Menon et al. , 1993; McMurtry et ai , 1993; Gao and O’Brien, 
1993; Madnia and Givi, 1993; Frankel et al, 1993; Cook and Riley, 1994; 
Givi, 1994; Fureby and Lofstrom, 1994; Moller et al. , 1996; Branley and 
Jones, 1997; Cook et ai , 1997; Jimenez et ai , 1997; Mathey and Chol- 
let, 1997; Colucci et al. , 1998; DesJardin and Frankel, 1998; Jaberi and 
James, 1998; Reveillon and Vervisch, 1998; Vervisch and Poinsot, 1988). 
Amongst these, recently Colucci et al. (1998) developed a methodology, 
termed the ‘‘filtered density function" (FDF). The fundamental property 
of the FDF is to account for the effects of subgrid scale (SGS) scalar fluc- 
tuations in a probabilistic manner. This is similar to probability density 
function (PDF) methods which have proven to be very useful in Reynolds 
averaging procedures (Libby and Williams, 1980; Libby and Williams, 1994; 
O’Brien, 1980; Pope, 1985; Dopazo, 1994). Colucci et al. (1998) developed 
a transport equation for the FDF in constant density flows in which the 
effects of unresolved convection and subgrid mixing are modeled similarly 
to those in “conventional” LES, and Reynolds averaging procedures. This 
transport equation was solved numerically by a Lagrangian Monte Carlo 
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procedure and the results were compared with those obtained by direct 
numerical simulation (DNS) and by a conventional finite difference LES in 
which the effects of SGS scalar fluctuations are ignored (LES-FD). It was 
shown that in non-reacting flows, the first two SGS moments of the FDF, 
as obtained by the Monte Carlo solution, are close to those obtained by 
LES-FD. The advantage of the FDF was demonstrated in reacting flows in 
which its results were shown to deviate significantly from those obtained by 
LES-FD but compare favorably with DNS data. The encouraging results 
generated bv FDF warrant its extension and application to more complex 
flows. Further assessment of its predictive capability is also in order. The 
primary objective in this work is to extend the FDF methodology for LES 
of three-dimensional (3D) turbulent reacting jet flows. The FDF deals only 
with scalar quantities; the hydrodynamic field is obtained via conventional 
LES. 


2. Formulation 


We consider constant density turbulent reacting jet flows involving N s 
species. The primary transport variables are the velocity vector u t (x, t), 
( i = 1,2,3), the fluid pressure p(x, £), and the species’ mass fractions 
<t>a Gl 0 {a — 1,2, N s ). These variables are governed by the conservation 
equations: 



(i) 


du , duiu, 

+ 3 

dp 

dr tJ 

dt dxj 

dXi 

dx 3 ' 


dJ 7 a 


dt dxj 

dxj 

+ 


(2) 


(3) 

where uj a is the chemical source term. Assuming a Newtonian fluid and 
Fickian diffusion, 


Tij = v 


( du x diij \ 

dxj dx{ J ’ 


r d<t>a 

dxj ’ 


(4) 


where v is the kinematic viscosity, T = ^ is the molecular diffusivity and 
Sc is the molecular Schmidt Number. Large eddy simulation involves the 
use of the spatial filtering operation (Aldama, 1990; Moin, 1991) 


(<!>(& t))i 



h s (x — x f )<t>(x f , t ), dx [ 


(5) 


where h $ (x) denotes the filter function of width A//, and (<f>(x,t))L rep- 
resents the filtered value of the transport variable <£(x, t). We consider 
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spatially & temporally invariant, localized and “positive* filter functions 
(Vreman et a/., 1994). The application of the filtering operation to the 
transport equations yields: 


d(uj) L d(uj) L (uj) L 
dt dxj 


9 (p)l . d{T {j ) L 

dx{ dxj 


d(4> a h , d( Uj ) L (4> a )L _ d(Jf) L dMf x 

~dT~ + dTj " + {Ua)L (8) 

where T tJ — { U{Uj)i - (u 1 )l(u j )l and Mf = {u 3 (j> a )i - denote 

the SGS stress and the SGS mass flux, respectively. 

The closure problem in LES of non-reacting flows is essentially one 
of representing the unresolved terms T x j and Mf. In reacting flows, the 
problem is compounded by the presence of the chemical source term 
for which an additional model is required. For closure of the hydrodynamic 
SGS stresses, the gradient-diffusion approximation is invoked: 

Tij ~ (ShMTi* = —2u t (Sij) L ( 9 ) 

where (S zj )l is the resolved strain rate tensor and v t is the SGS viscosity 

modeled via the modified kinetic energy model (MKEV) (Colucci et a/., 
1998): 

v t = CkA H ^\(u^) L (u’)i - «u*) L )/ / «M*) L ) z/ |, (10) 

where u * = u t — U x and U t is a reference velocity in the X{ direction. The 
subscript V denotes the filter at the secondary level which has a char- 
acteristic size (denoted by A H *) larger than that of grid level filter. The 
gradient-diffusion approximation is also used for closure of the SGS mass 
fluxes (Eidson, 1985): 

M ° = _ Tt 8ML (11) 


where T t = v t /Sct, and Set is the SGS Schmidt number and is assumed 
constant. 

The filtered density function (FDF) is utilized to represent the scalars 
in a probabilistic manner. For the scalar array (/>{xj) = [0i, • • -<f>N s ]i 

the FDF, denoted by /^, is defined as (Pope, 1990): 

= J e h 3 (x' - x)dx', (12) 


<£{*, 0] = % - £(*, 0] = if - 4>a(x,t)], 
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where S denotes the delta function and ip denotes the composition domain 
counterpart of the scalar vector <p. The term e[<p - ip(x,t)\ is the “fine- 
grained” density (Lundgren, 1967; O’Brien, 1980), and Eq. (12) states that 
the FDF is the spatially filtered, fine-grained density. Thus, Jl gives the 
density in the composition space around x, weighted by the filter h s . With 
a positive definite filter (Vreman et al ., 1994), fi has all the properties of 
the PDF. For further development, it is useful to define the “conditional 
filtered value” of the variable Q(x,t) by 

._ J-^Q(x',t)e \ip,<p(x', t)l h s {x' - x)dx' 

(Q(x, t)\rP) L = v -ff- -1 (14) 

_ fL{±,x,t) 

where (a\(3)i denotes the filtered value of a conditioned on (3. Equation 
(14) implies 


(*) For Q{x, t) = c, ( Q{x , t) \f) L = c 

(ii) For Q{x,t) =Q{<p(x,t)), (Q{x,t)\tp) L = Q(ip) (15) 

/ +oo 

(Q(x, t)\xp) L f L (ip ; x , t)drp 

-CO 

where c is a constant, and Q{<p{x, t)) = Q(x, t) denotes the case where 
the variable Q can be completely described by the compositional variable 
<p{x_,t). These properties, in conjunction with the FDF, facilitate the cal- 
culation of the moments involving the scalar variables via integration over 
composition space, 


/ +o° ^ 

Q(^)/t(^;x,t)d^. (16) 

-CO 

The FDF transport equation is obtained by taking the time derivative of 
Eq. (12) and making use of Eq. (3): 


dt dxj 


+ 


- ( uj)L]fL 


dx } 





d[2 a (±)f L ] 

dw a 


(17) 


This is an exact transport equation for the FDF. The last term on the RHS 
is due to chemical reaction and is in a closed form. The second term on 
the left hand side represents the filtered convection of the FDF in physical 
space and is also closed (provided («,•)/_, is known). The unclosed terms are 
the first two terms on the RHS which represents the transport of the FDF 
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via SGS convection and the effects of diffusion in composition space. The 
SGS convective flux is modeled via the gradient-diffusion approximation, 

[(u : \±)l - (uj)l]Il = -r (18) 

The closure for the conditional SGS diffusion is based on the linear mean 
square estimation (LMSE) model (O’Brien, 1980), which is also known 
as the interaction by exchange with the mean, or the IEM model (Borghi, 
1988). Implementation of this model together with Eq. (18) yields the mod- 
eled FDF transport equation: 


djL d[(uj) L f L ] 
dt dxi 


+ 


d_ 

dx{ 


(r + r,) 


dh 

dxi 


9 rr> \_NC .1 d lZ'a(±)fL] 

\4>a)L)fL] o , * (19) 




dy Q 


In the second term on the RHS , is the frequency of scalar mixing 
within the subgrid and is modeled via Q m = Cn(r + I\)/A^. This equation 
may be integrated to obtain transport equations for the SGS moments. For 
example, the first moment, (< P q )l > or the filtered mean is governed by: 


d{<t>g)L d(u 3 ) L {<t>a)L = 0 , d(<j) a ) L 

dt dxj dxj 1 dxj 


+ (^a)L? 


( 20 ) 


3. Numerical Formulation 

The numerical solution of the hydrodynamic and the scalar fields involves 
a two step explicit procedure. The first involves the advancement of the 
hydrodynamic variables and is accomplished via a compact finite differ- 
ence scheme (Kennedy and Carpenter, 1994). The second involves the time 
advancement of the FDF for which a Lagrangian Monte Carlo procedure 
is used. This procedure is based on the idea of “equivalent systems” by 
considering the random process A\(f), 

dXi(t) = Di(X(t),t)dt + E 1/2 (X(t), t)dW x [t), (21) 

where Di{X_,t) is the drift vector, E(X.,t) is the diffusion coefficient and 
Wi represents the Wiener-Levy process (Karlin and Taylor, 1981). With 
the equivalence: 


E = 2 (r + r t ), 


D t = (u,)l + 


<9(r + r f ) 

dx t 


( 22 ) 
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the Fokker Planck equation corresponding to stochastic differential equa- 
tion (21) becomes equivalent to the spatial transport equation of the mod- 
eled FDF equation (19). 

In the the numerical solution, the FDF is represented with a set of 
scalars t) assigned on the particles throughout the flow- field. 

The location of the notional particles are given by and Eq. (21) is 

integrated via the Euler-Maruyamma scheme: 

*i”Wi) = */ n) (4) + D i l n) (t k )At+ (23) 

where D\ n) (t k ) = D,(2C (n) (4), t), E^(t k ) = E(X (n \t k ),t) and £ t (n) is a 
random variable with the standard Gaussian PDF. This schemes preserves 
the Markovian character of the diffusion process (Gardiner, 1990) and fa- 
cilitates affordable computations. The coefficients D{ and E require the 
knowledge of the filtered mean velocity and the diffusivity (molecular and 
SGS). These are provided by the solution of Eqs. (6)-(7) by a finite differ- 
ence procedure and then is interpolated to the particle location. 

The scalar composition of each particle changes due to the effects of 
chemical reaction, and mixing (SGS and molecular). Both mechanisms are 
implemented deterministically and the scalars evolve according to 

^ = -n m (<PZ - ( 4>*)l ) + Won (24) 

where </>+ denotes the scalar value of a particle. 

4, Results 

Both FDF and LES-FD are employed for simulations of 3D turbulent round 
jets under both non-reacting and reacting conditions similar to those con- 
sidered in the experiments of Shea (1977). In the nonreactive case, the 
configuration consists of a jet of ozone (O3) diluted in nitrogen (N 2 ) issu- 
ing into a coflowing stream of N 2 - In the reacting flow, the surrounding fluid 
consists of nitric oxide (NO) diluted in N 2 . The chemistry is modeled via 
the one-step reaction of O3 + NO —>02 + N0 2 . The ratio of the reactants' 
concentration to that of the carrier gas is of order O(10~ 4 ). With such di- 
lute reactants, the effects of reaction exothermicity can be neglected (Shea, 
1977). In reacting flow simulations via LES-FD, the SGS scalar correlations 
are neglected. 

The streamwise velocity at the inflow boundary is initialized with an 
approximate top-hat radial distribution. The initial velocity is U Q in the jet, 
and Ucc, in the co-flow, with a velocity ratio of U 0 /Uoo = 4. The Reynolds 
number based on the jet diameter ( D ) and the inner jet velocity is Rep = 
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4, 000. The space coordinates are x — [x, y, z], where x is the streamwise di- 
rection, and y & z are the radial/cross-stream directions. A mesh consisting 
of 101 x 61 x 61 grid points in the x, y, z directions, respectively, is used to 
cover a domain of size 8 D x 4 D x 4D. The ratio of the secondary filter size 
to the grid filter size is Ah*/ A// = 3. The values of the other parameters 
are: Sc = 1, Set = 0.7, Ck — 0.045, Cn = 2. No attempt was made to 
find the optimum, or the “dynamically” determined (Germano et a/., 1991; 
Germano, 1996) values of the model constants. 

The simulation results are statistically analyzed via time averaging over 
16,000 samples. In the FDF simulations, the filtered values of the scalar 
quantities are determined by performing local averaging. The volume from 
which an ensemble of particles is constructed is A'g. By increasing A^;, the 
number of particles in the ensemble increases. This improves the statistical 
accuracy but increases dispersion. First, LES results of the non-reacting 
jet flow are considered in which the FDF simulations are conducted with 
Ae = 2A//. This size facilitates the use of fewer particles while still retain- 
ing a large enough sample for reliable statistics. The instantaneous density 
of the number of the Monte Carlo particles is presented in Fig. 1. This 
figure provides a visual demonstration of the basic methodology and the 
flow structure, as captured by the FDF. To establish the consistency of the 
FDF, its results are compared with those of LES-FD. Shown in Fig. 2 are 
the contour plots of the filtered ozone mass fraction at planes normal to the 
streamwise coordinate. The results via FDF are very similar to those ob- 
tained by LES-FD; the latter contain slight numerical oscillations which are 
not present in the Lagrangian Monte Carlo simulations. The comparison 
between the filtered values as predicted by FDF and LES-FD is quantified 
by performing a linear regression analysis of data at all the points. This 
analysis yields a correlation coefficient of 0.99 between the two sets of re- 
sults which indicates a very good agreement between the LES-FD and the 
FDF in predicting the filtered mean values. 

The radial d istrib utio ns of t h e time- averaged, filtered, normalized ozone 
mass fractions (Yo 3 ) L = (Yo 3 ) l/ (Yo 3 ) l( x = y = z = 0) are shown in Fig. 
3. In the non-reacting case, expectedly, the FDF results agree well with 
those via LES-FD. Both simulations predict a similar rate of decay for 
the centerline values of the mass fraction as the flow evolves. In the react- 
ing case, however, there is a significant difference between the results of 
the two simulations. It is noted that LES-FD predicts a much larger rate 
of reactant conversion in comparison with FDF. This difference is due to 
the neglect of the SGS scalar fluctuations in the LES-FD. This trend was 
observed in all the cases considered and is consistent with that observed 
in Reynolds-averaged simulations (Bilger, 1980). An attempt was made to 
compare the results with experimental data of Shea (1977). But there are 
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not sufficient detailed data reported in regard to the initial conditions in 
this experiments. Also, because of numerical concerns some of the simula- 
tion parameters are different from those considered experimentally. Work 
is in progress to generate DNS data for 3D, turbulent reacting jet flows 
with simple chemistry schemes of the type considered here. Such data are 
needed for further assessment of the methodology before it is implemented 
for simulations of more complex reacting flows. 

Although the FDF methodology is presented here for isothermal, con- 
stant density, reacting flows with a simple kinetics scheme, the extension to 
variable density flows, with exothermic reactions imposes no serious math- 
ematical difficulties (Jaberi et a/., 1999). For LES of variable density flows, 
it is convenient to use the filtered mass density function (FMDF), denoted 
by Fl, defined as 

F L {xjr,x,i) = J p(x',t)e h 8 (x f - x)dx', (25) 

where p is the fluid density. The integral property of the FMDF is such 
that 

/ +oo r+oo 

FL(}hx,t)dip= I p(x', t)h s {F - x)dx' = (p{x,t)) L . (26) 

-OO J — OO 

Jaberi et ai (1999) developed a transport equation for the FMDF and 
applied it for LES of several reacting flows. All the results as compared 
with DNS and laboratory data show significant advantages over LES-FD. 
With inclusion of efficient numerical integration routines for the treatment 
of complex chemistry mechanisms (Pope, 1997), it is conceivable that LES 
of reactive flows with realistic chemical kinetics may be conducted for en- 
gineering applications in the near future. In this regard, the scalar FDF 
methodology is attractive in that the present Monte Carlo solver can be 
used directly in available CFD codes. Similar to PDF methods, the closure 
problems associated with the FDF (and FMDF) are the correlations in- 
volving the velocity field (such as SGS stresses and mass fluxes). This may 
be overcome by considering the joint velocity-scalar FDF (FMDF) simi- 
lar to that in PDF methods (Pope, 1994b). This issue is currently under 
investigation. 

The computational requirement for FDF simulations with 2 x 10 6 parti- 
cles is about 2.5 times that of LES-FD. This overhead appears tolerable in 
view of the attractiveness of the methodology. Also, the computational re- 
quirements for FDF is significantly less than that of DNS. But the range of 
flow parameters (such as the Reynolds and the Damkohler numbers) that 
can be considered by FDF is significantly larger than can be treated by 
DNS, and the results are more accurate that those by LES-FD. Colucci et 
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al. (1998) and Jaberi et ai (1999) report a comparison of the computational 
requirements of LES-FD, FDF and DNS for several flow configurations. 
This comparison could be made only in flows for which DNS was possible, 
i.e. low Damkohler and Reynolds numbers. At higher values of these param- 
eters, the computational cost associated with DNS would be exceedingly 
higher than that of FDF. Thus for practical flows for which DNS is cur- 
rently impossible, the FDF would be a good alternative. Several means of 
reducing the PDF’s computational requirements are possible and should be 
considered. These could be useful in future applications in complex flows. 
The FDF method will benefit from ongoing and future improvements in 
PDF and other LES schemes (Pope, 1994a; Subramaniam and Pope, 1997; 
Pierce and Moin, 1998) from both modeling and computational standpoints. 
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Figure 1. Monte Carlo particle number density. 
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Figure 2. Instantaneous filtered mean ozone mass fraction contours at streamwise 
planes: (a) LES-FD, x/D = 2.5; (b) FDF, x/D = 2.5; (c) LES-FD, x/D = 7.5;(d) FDF, 
x/D = 7.5. 
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Figure 3. Time-averaged filtered mean ozone mass fraction. 
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Please Note: The materials provided in this Appendix (IV) are very preliminary and work is still 
in progress on this portion of our activities. 
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Velocity Filtered Density Function for Large Eddy Sim- 
ulation of Turbulent Flows 

I Introduction 

Modeling of the Sub-Grid Scale (SGS) correlations in large eddy simulations (LES) of turbulent 
flows is continuing to be an active area of research in fluid dynamics . 1-11 The most prominent 
model has been the Smagorinsky eddy viscosity closure 12 which relates the unknown subgrid scale 
Reynolds stresses to the local large scale rate of flow strain . 13 This viscosity is aimed to provide 
to a zeroth order approximation the role of mimicking the dissipative behavior of the unresolved 
small scales. The extensions to “dynamic” and “mixed” models 14-19 have shown some improve- 
ments. This is particularly the case in transitional flow simulations where the dynamic evolutions 
of the empirical model “constant” result in (somewhat) better predictions of the large scale flow fea- 
tures. More recent investigations and developments for more accurate SGS closures 20-27 imply the 
resolution of subgrid transport equations. The motivation in these cases is the desire to reduce the 
dissipative effect found in fixed length scale approaches and to predict the now famous problem of 
“back-scatter” by evaluating a proper velocity subgrid scale. Here again the “dynamic” approach is 
feasible to evaluate the required closure coefficients . 23 Two classes emerge in the subgrid velocity 
scale approach. The one equation type of approach investigated by Menon et al . 20 ' 21 solves for the 
subgrid kinetic energy, while the second one originally proposed by Deadorff et al . 25 and studied 
by Fureby et al P directly solve for the SGS transport equations. The latter model is theoretically 
more attractive as it has the ability of resolving anisotropy in the SGS. 

In a recent study Colucci et al. develop a new methodology for LES of turbulent reacting flows . 28 
In this new approach the definition of “Filtered Density Function ” 28 * 29 (FDF) and “Filtered Mass 
Density Function ” 30 (FMDF) are defined for the scalars and allow an exact representation of the 
chemical source terms appearing in the LES equations of turbulent reacting flows. While the ap- 
plicability of the scalar FDF for LES of chemically reacting turbulent flows has been successfully 
demonstrated by Colucci et al ., 28,30 the hydrodynamic part of the problem was limited to conven- 
tional Smagorinsky hydrodynamic closures 12 and of now well known limitations . 31,32 

The objective of the present work is to derive a higher order type of hydrodynamics closure through 
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the modeling of the joint- velocity FDF and to demonstrate its applicability by providing results 
based on its implementation for LES of turbulent flows. Assessment of the methodology is obtained 
at first through consistency simulations as explained later. Only the FDF of the joint-velocity vector 
is considered here; probabilistic treatment of the velocity-scalar fluctuations is postponed to future 
work but is the long term goal of the present investigation. 


II Formulation 


The primary transport variables for the mathematical description of incompressible (unit density) 
turbulent flows are the velocity vector u,(x. / ) (/' = 1 . 2. 3) and the pressure p(x. t ). The equations 
which govern the transport of these variables in space {x,) and time (/ ) are 


du , 
ot 


+ 



d 11,11, dp 
dx, d.r , 


+ 


dv,., 

dx , ' 


Assuming a Newtonian flow, the viscous stress tensor a :j is represented by 


( 1 ) 


IT, 


J 


V 


l. dx, dx, ) ■ 


where v is the fluid viscosity and is assumed to be constant. 

Large eddy simulation involves the use of the spatial filtering operation 33 


( 2 ) 


/ix. I 


-i: 


fix'. I )Q(x r . x)(lx . 


( 3 ) 


where Q denotes the filter function, (/(x. / ))/. represents the filtered value of the transport variable 
./lx./), and f = f — if) i denotes the fluctuations of / from the filtered value. We consider spa- 
tially & temporally invariant and localized filter functions, thus (7(x'.x) = <7(x' - x) with the 
properties, 33 ( x i = (!{— x), and (i(x)dx = 1. Moreover, we only consider “positive” filter 
functions as defined by Vreman etal . 34 for which all the moments x n, (i(.r)dx exist form > 0. 
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The application of the filtering operation to the instantaneous transport equations yields 


dxi 

djUjh d(u t )l.{' u :i)l. _ ^JPh tija,,),, _ rir L (i/ ; . Uj) 
tit ti.r, dxj dr, dr, 


( 4 ) 


where ti (v,. v, ) = {u,u ,)i — {»/,)/.(»,)/, denotes the “generalized subgrid stresses” as defined by 
Germano . 14 


Ill Deterministic modeling 

In LES the closure problem is associated with ri(u t . it, ). 2 Here, this term is modeled with a prob- 
abilistic and deterministic closures. The former is based on the velocity filtered density function 
and is discussed in the next section. For the later, several existing closures are given and used for 
comparison with the VFDF. 

The most famous deterministic LES model is probably the Smagorinsky model 12 which assumes 
equilibrium between the energy production and dissipation rates in the small scales. The model 
reads , 12 


T[. ( Uj. Uj ) — —2 U> Sjj + —hb, r 

_ 1 / d(u;) L ti(Uj)i.\ 

' " " A 07 , ) 

v, = A l S. 


( 5 ) 


r„i is a constant of order 0 . 01 , . S' = yJS~S~ and A/, is the characteristic length of the LES filter. 
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A more appropriate closure was proposed by Menon et al. Z0 ’ 21 in which the SGS are given by, 


d_ 

di 


[A-] + 




-2 Vt S tJ + -k 6ij. 


3 


= C»\ A/ V k. 


_0_ 

dxi 


.{ Uk)Lk ) = ~ 


t) 


d, 


t'k 


(v + I't) 


dk 

dxi. 


■ = c. 


kl 

a? 


Tl{u,.Uj) 


<Hu<) L 

d,fj 


( 6 ) 


and for which the introduction of k = 1/2 77, ( u a • « a ) as a subgrid velocity scale removes the equi- 
librium assumption allowing more realistic flow predictions. 

For our purpose we are interested in the transport equation of the SGS stresses, 14 

| b(«<- ";)] + ^ + p„ - ( 7 ) 


to allow resolution ofthe anisotropy present in the small scales. In this equation, T, jk = 7,(11, , Vj.u k )— 

Wjj] is the subgrid turbulent transport tensor where t l (u,.u r u k ) = (11,11,1^.) L -{u i ),T,\u l . u k )~ 
{"./}/ t,( u ,. u k ) - ( u k ) L TiXu< . ti ,) - The other terms are, the subgrid pressure- 

velocity scrambling tensor, 1 1 , = 77, ( a,. 


) + t l ( u j. j^), the subgrid production rate tensor, 




, and the subgrid dissipation rate tensor, 5 , , 


2 /'T L l 


"I U j 

(iTt ’ cfTi 


At this closure level a deterministic approach requires models for \ u,. a u k ), II, and 


Consistent with the models used in Reynolds averaged (RAS) calculations; the subgrid velocity- 
pressure scrambling tensor and anisotropic part of the subgrid dissipation rate tensor are combined 
and modeled via a Rotta type closure. 35 The resulting model is, 25,27 


■> 



■) 

Uj) ~ 7 k b,, . 


( 8 ) 


where ~ = j- is the subgrid mixing frequency, k = \t l (u,.u ; ) is the subgrid kinetic energy, and 
f = if,, is the subgrid kinetic energy dissipation rate. The dissipation rate is related to the charac- 
teristic length of the filter, A k , and the subgrid kinetic energy according to the same expression as 
in Eq. (6), 25,27 

The third order term, r, ( ti,. u , . a k ) , needs also to be modeled in a deterministic approach which is 
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( 9 ) 


done using, 25,27 


tl. r «*•) = - ^ 'i/2 -^L Vk 


_d_ 

dxi 


[ r L(w»- “j)j 


All the coefficients (\,C,,C U . have to be provided externally or may be calculated via dynamic 

methods. 14-18 Note also that more elaborate closures similar to those used in RAS could be used 
but they are beyond the scope of the present work. 


IV The velocity filtered density function (VFDF) 


The key point in this formulation is to consider the velocity vector U(x. 1 ) in a probabilistic manner. 
For that, we define the “ velocity filtered density function” (VFDF), denoted by P L , as 


P L { V: x.t) = J + X g [V. U(x'. /)] G’(x' - xWx'. 


g[V.U(x.l)} = S[V-U(x.1)} = Ipn; - M i( x -0] 

?=1 


( 10 ) 


where d denotes the delta function and V is the velocity state vector. The term £> [ \ 7 : U ( x. I )] is 
the “fine-grained” density, 36,37 and Eq. (10) implies that the VFDF is the spatially filtered value of 
the fine-grained density. Thus, / J y gives the one point, one time density distribution in the velocity 
space of the fluid state weighted by the filter G. With the condition of a positive filter kernel, 34 Pi. 
has all the properties of the PDF. 37 

For further developments, it is useful to define the “conditional filtered value” of the variable Q(x.t) 
by 


(Q(x. t)\V)i, = 


£2 g(x , ./) g [V.l7(x / ./)]G(x / -x)rfx / . 
Pi,( V: x. / ) 


( 11 ) 
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where (o 1 3) l denotes the filtered value of o conditioned on :3. Equation (11) implies 


(') For Q{x.l) = c. (Q{xJ)\V) L = c. 

( " ) For Q(x. /) = Q(U(x.1)). (Q(x. f)| V) L = Q{ V). (12) 

(Hi) Integral property : (Q(x.t)) L = f + (Q(x.1)\V) L P L {V:x.t)dV. 

J — (X 

where c is a constant, and Q(U(x.t)) = Q(x.t ) denotes the case where the variable Q can be 
completely described by the variable U(x. t ). Note that for simplicity the following abbreviation 
is used: ( .4 1 V 7 ) / = (A\U(x. t) = V)l- From these properties it follows that the filtered value of 
any function of the velocity variable is obtained by integration over the velocity space 


Q(x./)) l = J + " Q(V)P L (V:x.t )dV . 


(13) 


To develop a transport equation for the VFDF, the time-derivative of Eq. (10) is considered 

dPi(V-.x.i) /•<* dudx'.t) dg[V.U(x'.t)} 


0t 


L 


()i d\ 

<) r* diijix'. i 


C(x' — x) dx 


/ 

()\) ./_ 


(14) 


dt 


o [I 7 . U(x . /)] C(x' — x)dx'. 


This combined with Eq. (11) yields 

() Pi ( V : x. / ) <) 

Jfl = ~W, 

Substituting Eq. (1) into Eq. (15) yields 


(hi, t \ 

w 11 /, 


(15) 


()P,iV:x.1) _ () 

Jii “ ~di, 


'() 


u ha- 


th 


V 


n- 


' <h> 

Or, 




ft: 


n 


i j 


Pi.( V: x. / ) > . 


(16) 


in which the convective term can be represented in the form 


c) 

JJv 


()u, 


<n 


<).v k 


P,(V:x.l : 


u 


. ()P/ (V: x. /) 


().V k 


(17) 
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Substitution of Eq. (2) for the conditional diffusion and some algebraic manipulations yield, 




Pl(V:xJ) 

L 


d 2 P L (V:x.t) d 2 

// -|- 

dx k dx k dVjdVj 


/ ()u , da 

/ y 

\ dx k d.r k 


P L (V:x.t) . 


(18) 


The conditional filtered values of the pressure gradient and dissipation rate can be further decom- 
posed into resolved and deviatory components from the mean (filtered) values. It is useful to adopt 
the following decompositions, 


1 k Pi. = (u^lPl + [V* — (u k )l]Pl- 


~\v) Pl = ^ ± Pl + 

dx, / Oxi 


&iv\ -Mi 

dx, dx i 




dx k dx k 


dx k dx k 


du^^yX p (19) 

dx k dx k j L <>x k dx k \ 


so that Eq. (15) can be expressed as 


() uv hPll d(p) L dP L d 2 p L d(v,) L d{u,) L d 2 p L 
fx k dx, d\ , dx k dx k dx k ()x k d\ ,d\ , 


d \ \\d 


d\]d\ J 


duj du , 
/ 

dr k dx k 


V) - 


d(u,) L d(u,)r 
dx k dxi. 


Pl • 


where jjj = + (</*.) h denotes the filtered material derivative. 

Equation (20) is an exact transport equation for the joint- velocity FDF. The first term on the right 
hand side represents the deviatory /subgrid convection of the VFDFin physical space and is closed 
(provided that ( u k ) / is known). The second term corresponds to the convection in the velocity space 
due to the resolved pressure gradient. The third term is the diffusion of the VFDF in physical space 
due to the molecular effects, and the fourth term is the diffusion in velocity space due to the resolved 
dissipation rate. The unclosed terms are associated with the last two terms on the RHS of Eq. (20). 
These terms represent the convection in velocity space by the unresolved subgrid/deviatory pressure 
gradient and the diffusion in velocity space by the unresolved subgrid/deviatory dissipation rate. 

The subgrid pressure gradient and the subgrid dissipation rate are modeled via the generalized Langevin 
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model 38-40 


d 

d\) 



<Hph \ 

dr, ) 



d 2 

dVidVj 



\\ dr, dr, 



d(uj) L \ p 
' dr, c)r, ) 1 


[Gij(V 3 -(vj) L )PL] + -C 0 c 


d 2 Pi 

dv~m 


( 21 ) 


The advantage of the decomposition in Eq. (20) and the subsequent model in Eq. (21) is that they 
yield results “nearly” similar to those in “conventional” LES for the first two moments of the VFDF. 
To show this mathematically the moment equations are evaluated by integrating Eqs. (20 & 21) in 
velocity space. These moment equations are read as, 


djUj) L 

dr, 

d{u,)L d(u,) L (u,)i 


0. 


di 


dr, 


d(p)L + ^(vjh 


dxj 


dr,dr, 


dr L (ui. u,) 
dr, 


( 22 ) 


Uj)J + ~-l(u t ) L r L («,. »,)] = 


c) c) 

blK- «*) - >'T7— (/,)]] 

UZ'k 


+ P ik T l{ U r Ok) + (jj k T I,(o,. U,) — T/Jll,. U,) ■ ^ — — T]{Vj. U ,) — ' ' L 


dr, 


dr, 


+ Coshr 


(23) 


The advantage of the VFDF approach is seen in Eqs. (22,23) where the subgrid stresses appear in 
a closed form without the need of solving subgrid transport equations. Moreover, the third order 
quantities, [u,. u . . u, ), appears as a consequence of the VFDF model and do not need to be mod- 
eled separately as in a deterministic closure, i.e. Eq. (9). The accuracy of the subgrid stresses as 
obtained from the modeled VFDF transport equation for Eq. (20) need nonetheless to be validated. 

To make the second moment equation derived from the VFDF transport equations more similar to 
the one derived directly from the instantaneous equation (7), we define , as, 41 


(A, - 


i - 
•) 




) <\ 


.£ 1 , 

** .y 


(24) 
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With this, the resulting transport equation for the modeled VFDF becomes. 


DPl 

1)1 


3 



-[(Vi- -(«*) l)/ j l] + 

k 


d{p) L VPl 

().r, () \ 


— v 


Gij (V/ ~ ( u j)l) I J L ] + - Co 


d 2 P L 

dv,d\ 


d{tn)L 
dx k dx k 


d 1 Pi 

d\,d\j 


+ V 


d 2 P L 

().r k d.r k 


( 25 ) 


There are two constants in Eq. (25). While the first one, C e , is expected to be filter size dependent, 26 
the second one, C u , is equal to 2. 1 for very high Reynolds number flows according to Obukhov- 
Kolmogorov hypotheses. 42-46 For finite Reynolds number flows the value of C 0 can vary between 
0 and 6 based on theoretical works and RAS simulations. 40,47 " 49 


V Monte Carlo solution of the VFDF 

The solution of the VFDF transport equation (Eq. (25)) provides all the statistical information per- 
taining to the velocity vector U(x.t). This equation can be solved most effectively via the Monte 
Carlo schemes which can be utilized in both Eulerian 50 and Lagrangian 40,51 contexts. The advan- 
tages of Lagrangian numerical methods in reducing the amount of numerical diffusion are well- 
recognized. 52-57 The basis of the Lagrangian solution of the VFDF transport equation relies upon 
the principle of equivalent systems . 37,41 Two systems with different instantaneous behaviors may 
have identical statistics and satisfy the same VFDF transport equation. In the Lagrangian Monte 
Carlo procedure each of the particles obeys certain equations which govern its transport. These 
particles undergo motion in the physical space by convection due to the velocity vector and diffu- 
sion due to the molecular viscosity. The evolution of the velocity vector is governed by the cumu- 
lative effects of the local pressure gradient and velocity dissipation rate at the resolved and subgrid 
scales. The spatial and velocity diffusion of the particles are represented in a stochastic manner by 
the following system of stochastic differential equations (SDE) 37,58-60 

dX, (t) = /;,(X(/).V(/):/| (It + P(XU). V(7):/) d\\J(l). 

(1C, (I) = .W,(X(/).V(/):/) (II + L(X(/).V(/):/) JU7(/) (26) 

+ F„(X(/).V(/):/) 
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where A', is the Lagrangian position of a stochastic particle with velocity ! The coefficients I), 
and M, are known as the “drift” in spatial and velocity coordinates respectively. B and K are the 
“diffusion” coefficients for physical and velocity spaces respectively, and 1 1 and IT/ denote inde- 
pendent Wiener-Levy processes. 61 The tensor F,, represents the dependency between the velocity 
and physical spaces for this process. 

The corresponding Fokker-Planck equation for this Lagrangian process is. 


w + Jk lwv ' X)n - ~m m,v ’ X)r] + '2^ l (flJ|V ’ X)r 


d 2 


+ - 


2dv k dv k 
1 d 2 


(E 2 \V,X)p 


+ 


d 2 


dx k d \ ; 


[(Fa-B\V,X)p 


•>d VidVj 


l(F lk .F jk .\V,X)r). 


(27) 


where /* = /' ( V . X : / ) , is the PDF of the process, ( .4 ) is the “ expected value of A ”, and ( .4 1 V, X ) = 
{ .4 1 U ( / ) = V. X(1) = X) . The PDF of X , //, is obtained by integrating Eq. (27) over the velocity 
space, 58,62 




(28) 


Proper initial and boundary conditions for Eq. (28) ensures that /;( X . I ) = ( ’ xt (non zero). If this 
is satisfied one can divide Eq. (27) by /’ and derive the Eulerian transport equation. 60,63 

By comparison of the Eulerian transport equation corresponding to Eq. (26) with the modeled VFDF 
transport equation (25) one can determine a set of appropriate values for the coefficients. For ex- 
ample. 


.1/ = 


V (;>)/. 

()x, 


() 2 ( U ;) 1 

2l '~ — 1 — b (/,, (l j — {(/,)/,). I), : l 


().l k (). V k . 


B = \2u. F = F,j = \F>1> 


<){">) 
().v , 


(29) 


- p t ) L = {F,). 


is one convenient set of relationships but is not unique. With the equalities given in Eq. (29) /* 
is non zero as it should be, if the proper initial and boundary conditions are applied. Note that the 
diffusion coefficient, F = P2v is selected consistently with the one used by Einstein, 64 Wiener 65 
and Levy 66 to describe Brownian motions. Thus the chosen SDE’s which represent the transport of 
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the modeled VFDF are 


dXi(t) = r,(t)dt + yj'lv dW*(t). 

dUAt) = 


<Hp)l , n . d 2 (u,) L , ^ / \ x 

+ — 4 G,j (l .,(/) - 


+ 


dxkdx k 

d\V*(t). 

OX, J 


di + dcTeM'nn 


(30) 


This stochastic system reduces to the one developed for RAS calculations 60,67 when the filter oper- 
ation is assimilated to an ensemble averaging operation. 

It has been shown by Pope 37 that the stochastic system yields identical statistics as the fluid particles 
(to a first order approximation in time) if and only if the time increment, dt, is much smaller than 
the characteristic time of the large scale turbulent motions, r, and much larger than the Kolmogorov 
time scale r n : 


t„ < df <C r. 


(31) 


Under the previous constraints and to a first order approximation in time, the statistics obtained 
from the stochastic process evolve as: 


(d. V,) = (v,)i, dt. (dl'i) = 
(dX,d.\ t) = 2/' <S (/ dt. (df ,dt ,) 
(. dX-.dV ,) = - dt. 

()x t 


X{m) L 

— b 2.V 

0.V, ()-V k ()x k 


dt 


().r k . ().v k 


dt 4 C oc b,j df. 


(32) 


In the numerical implementation, the VFDF is represented by a set of Monte Carlo particles, each 
with a velocity vector / ) and a Lagrangian position vector X (n) . The simplest means 

of simulating Eq. (30) is via the Euler-Maruyamma approximation 68 

A 7 (/*. + i ) = A 7 ( 4 ) 4- ir;n,. )A / 4 /r (/*.)( A/ ) l/2 c;'(/*). 

r,"(/ A . +I ) = r;(/,i 4 .u"(/,)A/ 4 /;••(/. mA/) ,; - , ci/,) (33) 

+ /•;;;(/ a )(A/) i/2 c;(/ a .). 

where D"{t k .) = !) l (X { " i (l k ). U , " , ( (/* ): / ), /*"(/*.) = /y(X ( ">(/, ). . and £■' (t /,)•("(( 
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are samples from two independent Gaussian white noises at time / k . This formulation preserves the 
Markovian character of the diffusion processes 69-71 and facilitates affordable computations. Higher 
order numerical schemes for solving Eq. (30) are available, 68 but one must be very cautious in using 
them for LES. 28 Since the diffusion term in Eq. (26) strongly depends on the stochastic processes 
U(/) and X(/), the numerical scheme must preserve the Ito-Gikhman 72,73 nature of the process. 
Equation (33) exhibits this property. 

The LES of the first hydrodynamic moments is conducted with the “compact parameter” finite dif- 
ference scheme of Carpenter. 74 This is a variant of the McCormack 75 scheme in which a fourth 
order compact differences are used to approximate the spatial derivatives, and a second order sym- 
metric predictor-corrector sequence is employed for time discretization. The computational scheme 
is based on a hyperbolic solver which considers a fully compressible flow. Here, the simulations 
are conducted with a low Mach number (M % 0 . 3 ) to minimize compressibility effects. The proce- 
dure involved in the finite difference discretization is dependent of the Monte Carlo solver through 
the SGS. All the finite difference operations are conducted on fixed and equally sized grid points. 
Thus, all the filtered values of the hydrodynamic variables are determined on these grid points. The 
transfer of information from these points to the location of the Monte Carlo particles is conducted 
via interpolation. Second-order (bilinear) interpolation scheme is considered, as no significant dif- 
ference in statistics were observed when higher orders were used. 29 

The SGS necessary to the finite difference LES solver (or higher order moments of the VFDF) at 
a given point are estimated by consideration of particles within some volume centered at the point 
of interest. Effectively, this finite volume constitutes an “ensemble domain” characterized by the 
length scale A / in which the VFDF is represented discretely by stochastic particles. This is neces- 
sary as, with probability one, no particle will coincide with the point as considered. 40 Here, a box 
of size A / is used to construct the ensemble mean variances and covariances of the velocity vec- 
tor from which the ensemble mean values are subtracted to yield the SGS at the finite difference 
nodes. These values are used in the finite difference LES solver of Eq. (4). The subgrid kinetic 
energy dissipation rate and the subgrid mixing frequency are also obtained from the SGS. From the 
numerical standpoint, determination of the size of the ensemble domain is an important issue as it 
determines the time evolution of the LES solver through the values of the SGS. Ideally, it is desired 
to obtain the statistics from the Monte Carlo solution when the size of sample domain is infinitely 
small (i.e. A/ — > (I) and the number of particles within this domain is infinitely large. With a finite 
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number of particles, if Af is too small there may not be enough particles to construct the statistics. 
A larger ensemble domain decreases the statistical errors, but may increase the dispersion errors 
which manifest themselves in “artificially diffused” statistical results. This compromise between 
the statistical accuracy and dispersive accuracy as pertaining to Lagrangian Monte Carlo schemes 
implies that the optimum magnitude of Ae cannot, in general, be specified a priori . 31 This does not 
diminish the capability of the procedure, but exemplifies the importance of the parameters govern- 
ing the statistics. A better understanding of the sample size impact is obtained through consistency 
simulations as illustrated in the results section. 


VI Results 

VI.1 Flows simulated 

In this section results are presented to demonstrate the effectiveness of the VFDF method. Spatially 
developing jet configurations are considered for the LES simulations. 21) planar jet simulations are 
used for consistency assessment of the previously described methodology and for comparisons with 
existing closures. ID round jet simulations allow validation of the new approach and pre-existing 
closures via experimental data. 

All of these flows are dominated by large scale coherent structures. The formation of these struc- 
tures are expedited by imposing low amplitude perturbations at the inflow boundary. In the figures 
presented below, x. y correspond to the streamwise and cross-stream directions, respectively. In 
3D, r denotes the spanwise direction. Finally r = v //r + ~ 2 is the radial direction. The size of the 
domain in the 2D jet flowis 0 < x < 1 \1). — -iA D < y < I ). The ratio of the inlet jet velocity 
to that of coflowing stream is kept fixed at OA for the consistency analysis and 0.2 for comparative 
results between closures. For the :j D jet the domain consists of a rectangular box of dimensions 
0 < x < I I/). < y < 3M). 

All the flow configurations are simulated via LES. The procedure in LES is based on the Monte 
Carlo solution of the modeled VFDF transport equation (Eq. (25)) for the velocity vector augmented 
by the finite difference solution of the filtered hydrodynamic modeled equations (Eqs. (4)). In the 
presentation below, these results are identified by the abbreviation VFDF. In addition, another LES 
is conducted in which the modeled transport equations for the filtered velocity and the generalized 
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subgrid stresses are simulated with the finite difference scheme. In these simulations, the hydro- 
dynamic solver and the models for the SGS are consistent with those employed in the VFDF (i.e. 
Eqs. (7), (8) and imposing Eq. ( 24 )). The effects of the turbulent subgrid transport, T[ (u,. u r u k ), 
are obtained from the VFDF in the case of the consistency simulations. The results based on this 
procedure are referred to as LES-FD. No attempt is made here to determine the magnitudes of the 
constants appearing in these models in a dynamic manner. 14 However, different values are consid- 
ered for ( o in order to compare the predictions obtained from the VFDF with other deterministic 
closures. 

VI.2 Numerical specifications 

The primary parameter is the flow Reynolds number (Rt ). All finite difference simulations are con- 
ducted on equally-spaced square grids (A x = A y = A.:(for 3D ) = A). Since the size of the com- 
putational domain is fixed, the number (and the size) of the grids depends on the type of simulation 
being conducted. 

The VFDF and LES-FD runs are conducted on grids coarser than those in DNS. Unless otherwise 
specified, the LES resolutions in the consistency simulations of the planar jet are 201 x 101 and 
181 x 91, with R ( = 4.000 (based on the inlet jet diameter). The planar jet configurations used 
for comparisons of the various closures use 101 x 81 points with R< = 10.000. The :W grid is 
composed of 1 0 1 x 8 1 x 8 1 grid points for a Reynolds numbers of R< = 10. 000 based on the inner 
stream velocity and jet diameter. A low speed coflow corresponding to a 0.2 ratio with respect to 
the inner flow is maintained in order to stabilize the solver. 

When required (inlet conditions, DNS) a top-hat filter function 33 of the form below is used 

■V/. _ 

(7(x' - x) = JJ (',{■>', - .r, ) 

/ = 1 
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in which .\ r , denotes the number of dimensions, and A f , = 2 A. No attempt is made to investigate 
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the sensitivity of the results to the filter function 34 or the size of the filter. 76 

In VFDF, the Monte Carlo particles should initially be distributed at / = 0 throughout the domain. 
For the jet simulations the particles are supplied only in the region - 1 .75 D < y < 1.75 D in 'll). 
This procedure seems to be sufficient to yield accurate results, at least for the velocity ratios under 
study. The particle density is monitored at all times to ensure an approximately uniform particle 
density through the domain of interest. 

In the spatial jets, new particles are introduced through the inlet boundary at a rate proportional 
to the local flow velocity and with a compositional makeup dependent on the y, (and r in 3D) co- 
ordinate and yielding identical statistics regardless of the ensemble domain size A £ . The density 
of the Monte Carlo particles is determined by the initial number of particles per grid cell ( V PC) 
of dimension A x A( x A). The magnitude of X PG is varied to evaluate its affect on statistical 
convergence of the Monte Carlo results. This assessment is demonstrated in 2D simulations of the 
spatially developing planar jet. The simulations of 3D spatial jet are based on XP G = 40. Thesize 
of the “ensemble domain” in the VFDF simulations is also varied in order to quantify its influence 
on the statistical convergence. The following sizes are used, A/, = 2 A. A. A/2 in the consistency 
simulations and X E = X otherwise. The number of sample particles used to construct the VFDF 
statistics is thus controlled by the values of X PG and A^. 

An additional parameter which influences the numerical accuracy is the magnitude of the incremen- 
tal time step. The stability criterion for the finite difference scheme requires GPL < 1 / vA 74 and 
is more stringent than the criterion for the Fourier number. The effect of the time increment on the 
accuracy of the Euler-Maruyamma scheme is not investigated here. The A / value (CFL numbers) 
equal to the value from Colucci et al , . 28 is adopted. 

The simulated results are analyzed both “instantaneously” and “statistically.” In the former, the 
instantaneous contours (snap-shots) of the vorticity and scatter plots of the redundant variables are 
considered. In the latter, the “Reynolds-averaged” statistics constructed from the instantaneous data 
are considered. In these spatially developing flows this averaging procedure is conducted via sam- 
pling in time. All Reynolds averaged results are denoted by an overbar. 
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VI3 Consistency of VFDF and convergence of the Monte Carlo procedure 


The objective of this subsection is to demonstrate the consistency of the VFDF formulation and 
the convergence of the Monte Carlo simulations. For this purpose, the LES results via VFDF and 
LES-FD are compared against each other in a planar-jet configuration under different conditions. 
Since the accuracy of the finite difference procedure is well-established (at least for the first order 
filtered quantities), this comparative analysis provides a good means of assessing the performance 
of the Monte Carlo solution of the VFDF and allows identification of the governing parameters. For 
simplicity and a clear understanding of the methodology, the model is simplified by neglecting the 
spatial diffusion of the VFDF due to the molecular term. This assumption results in the neglect of 
the molecular transport in the SGS equations (Eq. (23)). The Generalized Langevin model in turns 
models the expression, 
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This simplification is equivalent to an assumed high Reynolds number at the subgrid level since only 
Eqs.(22) are fully recovered. This approximation reduces the degrees of freedom of the stochastic 
process. The new equivalent stochastic system therefore reads, 
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and has statistics evolving (to a first order in time) as, 
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The third order term in Eq. (7), r L (a,. u r o k ), is obtained from the VFDF and used into LES-FD. 37 
This ensure a fully consistent approach between LES-FD and the VFDF. The model’s coefficients 
are taken to be ( \, = 1 48 and (\ = \ . 77,78 Eq. (36) model could be used for LES in particular cases 
but the reader has to recognize that the proper behavior of the VFDF is not fully recovered. 79 

In Fig. 1, results are presented of the LES of the spatially developing planar-jet. Shown in the figure 
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is the instantaneous vorticity contour plots via (a) VFDF and (b) LES-FD. This figure provides a 
simple visual demonstration of the consistency of the VFDF as the results via the particle method 
are in agreement with those obtained by LES-FD. In fact, the Monte Carlo results are even more 
attractive due the Lagrangian nature of the solution procedure. While the LES-FD results suffer 
from slight over- and under-shoots, there are no such errors in the Monte Carlo scheme. The vi- 
sual agreement is confirmed by the scatter plots of the redundant variables (u) L and (r) L shown 
in Fig. 2(a) and (b). The correlation and regression coefficients for the first filtered moments were 
found to be rather insensitive to A E (±0.5%) as seen in Tables (1) and (2). 

Another rigorous means of assessing the VFDF results is via consideration of the Reynolds aver- 
aged results. Figure 3 shows such results in which the sensitivity of the VFDF predictions to sev- 
eral parameters is assessed. Figure 3(a) shows the comparison of VFDF and LES-FD results for 
the cross-stream variation of {u)j_ at a distance of 7.5 inlet diameters from the inflow boundary. 
The results are qualitatively identical at all y locations in the flow field. Several values of A e are 
taken under consideration for a unique number of particles per ensemble (X PC,' = 10.40.1 60). It 
is shown that the first filtered moment of the VFDF agrees very well with those obtained by LES- 
FD, even for large A f values. The differences between the VFDF and LES-FD results are more 
significant in Fig. 3(b) where the stream-wise variation of {</)/, is shown for several values of A / . 
This figure also indicates that the difference between VFDF and LES-FD predictions diminishes as 
Ay decreases. 

The differences observed in the first order filtered quantities are directly correlated to the contri- 
bution of the subgrid second order filtered moments. Significant differences observed at various 
instantaneous times between the subgrid moments from LES-FD and from the VFDF would result 
in distinct evolutions of the flow field in time. It is therefore necessary to have a better understand- 
ing of the effect of Ay on the subgrid quantities. Figure 4(a,b,c) shows scatter plots at a given time 
for (a) -y. ((/.(/), (b) ry ( r. r) and (c) -/.(</. r), for various A/. The corresponding correlation and re- 
gression coefficients are shown in Table (1). Convergence at the subgrid level is clearly obtained as 
Ay — * 0 . Note nonetheless the large difference between LES-FD and VFDFwhen Ay, = 2A. Time 
averaged quantities corroborate these last results as observed in Fig. 5(a,b,c) where (a) t l (u. v), (b) 
ry ( r. r ) and (c) -y ( i/ . r ) are shown. 

The other parameter which influences the accuracy of the Monte Carlo results is the number of 
Monte Carlo particles per grid cell (\ PC). Figures 6(a,b) and 7(a,b,c) show that {</)/. and JvjJ. 
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are totally insensitive to V PC for the range considered and that 77, (v t . u , ) values do not vary sig- 
nificantly for X PC > 40 with a given At. More important figures 6 and 7 illustrate the greater 
influence of the size of the ensemble domain than A PC. Consistency simulations with no third 
order subgrid terms in the LES-FD scheme were conducted and significant drops in the correlation 
coefficients were observed especially for At = 2 A (cases not shown). The Reynolds averaged 
subgrid stresses did not depict such drastic behaviors but differences could still be observed stress- 
ing the importance of the third order terms for any LES-FD/deterministic approach. This problem 
is avoided with the VFDF. Therefore in the general case even smaller X PC values can be used as 
long as At is properly chosen. This last statement is advanced based on the observation that the 
third order subgrid quantity v u*) is required for consistency simulations which is not the 
case for the general approach. 

The convergence of the Monte Carlo solution and the independence to X PC and At are demon- 
strated by these results (at least for the first filtered moments). The size of the ensemble average is 
nonetheless important as it influences the general behavior of the LES solver through the predic- 
tions of the SGS. It seems necessary to keep A t reasonably small in accordance with the theoretical 
point of view to estimate the subgrid quantities correctly. 

VI.4 Qualitative study of the VFDF: comparison with existing closures 

Deterministic closures were presented in Section II. The first one referred to as the “Smagorin- 
sky model”, Eq. (5), assumes equilibrium between the production and dissipation rates at the sub- 
grid scales yielding a zeroth order closure. The second closure, Eq. (6), referred to as the “A -eqn 
model”, 20 suppresses the previous equilibrium assumption by solving a modeled transport equation 
for k- = -jiX u k . u k ) which, to some extent, allows resolution of the present desequilibrium in the 
subgrid scales. Both of these models are nonetheless unable to resolve the anisotropy of SGS which 
is expected to increase in low resolution LES simulations of practical engineering flows. The last 
deterministic model corresponding to Eqs. (7,8,9) and referred to as the “SGS-eqns model” 25 - 27 the- 
oretically predicts anisotropy of the SGS as well as the desequilibrium process. The VFDF model, 
Eq. (30), has all the properties of the “SGS-eqns model”. These last two models are not fully con- 
sistent in this subsection. The third order quantities, -/ (u,. u being implied in the VFDF clo- 
sure have to be modeled with “SGS-eqns model”, Eq.(9). These term contributions were already 
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found to be critical in the previous subsection and are further studied for the second order closures. 
The coefficient Co found in the VFDF and “SGS-eqns” models is varied in order to evaluate its 
effect on the flow predictions. Only the Reynolds averaged quantities are considered here as they 
are illustrative of the instantaneous behavior of the studied dynamic systems. Our objective is to 
illustrate through numerical simulations the fundamental differences between the dynamic systems 
composed of the LES governing equations plus the LES model and quantify the closure coefficient 
significance. 

Figure 8(a) shows the values of ( it ) / along the center line as obtained from the different models. The 
Smagorinsky and A’-eqn models predict roughly the same slope of decay of the streamwise velocity. 
The higher order models on the other hand predict different behaviors depending on the value of 
Co- This observation allows us to stipulate on the importance of this coefficient in predicting the 
strength of the diffusive effect of the SGS scales. Also to notice from Fig. 8(a) are the different 
locations of the starting velocity decay from one model to another. The deterministic models predict 
the transition to occur at x = 7.5 D ± 0.5 D. The VFDF produces two distinct regimes. The first 
regime depicts a slowly decaying (v) L until r ~ 8/) where transition to a fast decay occurs. In the 
second regime the slopes are similar to the values obtained with the SGS-eqns closures. Shown in 
Fig. 8 (b) is the cross-stream variation along ,t = SI) for the streamwise velocity component non- 
dimensionalized by its value at the center line. Clearly the level of diffusion in the profiles decreases 
as one uses the Smagorinsky, the k -cqn model, the SGS-eqns or VFDF models. Note that results 
obtained without model are adjoined to the figures to illustrate the importance of the SGS in the 
present simulations. 

The decay rate of the center line velocity is related to the large structures behavior which is essen- 
tially governed by the subgrid production rate of kinetic energy, , . u L ) . This last quan- 

tity transfers the mechanical energy contained in the large structures to the internal energy and is 
dissipated by the viscous forces. The Reynolds averaged production rate and subgrid kinetic energy 
are shown in Fig. 9(a,b). The choice of LES model defines the mechanism of energy tranfer between 
scales. The Smagorinsky and / -eqn models predict roughly the same level of energy transfer while 
the other models predict more than twice the amount of subgrid kinetic energy production rate. Here 
again ( ' () is found to be an important parameter for the second order type of closures (SGS-eqns and 
VFDF). For approximately the same levels of subgrid turbulent kinetic energy, Fig. 9(b), the sub- 
gnd production rate differs quite much for various ( ‘ (J values. In Fig. 9(a,b) the subgrid quantities as 
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obtained from the VFDF models show two regimes as observed in Fig. 8(a). The subgrid produc- 
tion rate and turbulent kinetic energy first grow slowly at two different rates until x =s SI). Within 
this potential core energy is slowly removed from the large scales increasing the energy in the small 
scales. This transfer of energy results in a slowly decaying center line velocity component as ob- 
served in Fig. 8(a). At x % SD, k and its production rate suddenly increase to the levels obtained 
with the SGS-eqns model. This two regime behavior observed with the VFDF in Figs. 8(a)-9(a,b) 
can only find its source in the third order term as it is not predicted with the deterministic mod- 
els which ignore t l ( Uj. u k ) or supply a more or less ad hoc modelisation. Finally, differences 
observed between VFDF simulations can only result from the values of ( ' 0 as .V PC = 10 and 
A e — A are kept identical. 

Equation (23) demonstrates that C 0 influences only the anisotropy of the SGS. Hence, the obser- 
vations formulated from Fig. 9(a,b) are explained through the level of anisotropy in the SGS as 
predicted by the models. The third order term, tl(u;. u r u k ), can aslo be at the origin of these dif- 
ferences, but only to some limited extent. Anisotropy in the normal SGS are shown in Fig. 10(a,b) 
for the cross-stream direction along .?• = SD. The effect of C 0 is clearly illustrated at the subgrid 
level in these last figures. Because of the importance of the SGS in Eqs. (4) and in the evolution 
equation of k, Fig. 8, f 0 plays a crucial role in predicting the desired behavior of the LES dynamic 
system. 

The final part of this subsection considers the VFDF equations obtained in Eq. (30) and Eq. (36). 
These models are respectively referred to as “VFDF1” and “VFDF2” in the figures. The compar- 
ison is conducted for the same conditions as above and aim to study the differences between the 
two stochastic systems given in Eqs. (32) and Eqs. (37). The two models essentially differ in the 
presence or not of the resolved disspation rate in the time evolution of the incremental velocity cor- 
relations. Figure 10(a,b) validates the approximation made in Eq. (36) for the flow configuration 
studied. The resolved part of the dissipation rate is still found to have effects on the predictions and 
use of the original model, Eq. (30), is still advised as LES usually deals with M), transient flows of 
finite R< number in which resolved scale contributions are of critical importance. 

The comparative study of the behavior of various LES systems (i.e: LES solver plus SGS model) 
underscores fundamental differences in the dynamics of the systems inherited from the SGS mod- 
els. More specifically, the third order filtered correlations are found to be of importance for the 
second order closure systems. The C„ coefficient, necessary for the VFDF closure, seems critical 
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in estimating the anisotropy in the SGS as well as the proper decay rate and energy transfer between 
scales. The choice of coefficient in the VFDF model will therefore result in a specific dynamic evo- 
lution of the LES solver as observed in this subsection. As a consequence, gathering informations 
about Co for different flow configurations is necessary for validation of the VFDF approach. 
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Figures and Tables 
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Figure 1 : Instantaneous snapshot of the vorticity field: (a) VFDF, (b) LES-FD. 
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Table 1: Regression, r EF , and correlation, pi;r, coefficients as a function of Af for Fig. (4). Pla- 
naijet configuration: 201 x 101. 
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Table 2: Regression, />/., and correlation, p Ef ., coefficients as a function of A/ . Planaijet config- 
uration: 1-81 x 91. 
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Figure 2: Scatter Plots, /> / and pu are the regression and correlation coefficients, respectively - 

(a) {")/-> (b) ('•)/.. 
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Figure 4: (a) u scatter plot, (b)-/.(r. r) scatter plot and (c)-/.( u. r ) scatter plot. 
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Figure 7: Cross stream variations of (a) r L (//.//), (b) r L ( v . r ) and (c) t l ( u . v ) as a function of N PC 

{r = 5/X 9/7). 
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Figure 8: (a) Streamwise and (b) cross stream variations at r = SI) for ( a) / obtained with various 
LES models. 
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for Turbulent Reacting Flows 
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Explicit algebraic scalar-flux models that are valid for three-dimensional turbulent 
flows are derived from a hierarchy of second-order moment closures. The mathematical 
procedure is based on the Cayley - Hamilton theorem and is an extension of the scheme 
developed by Taulbee. Several closures for the pressure- scalar gradient correlations are 
considered and explicit algebraic relations are provided for the velocity - scalar correla- 
tions in both nonreacting and reacting flows . In the latter, the role of the Damkohler 
number is exhibited in isothermal turbulent flows with nonpremixed reactants. The rela- 
tionship between these closures and traditional models based on the linear gradient- 
diffusion approximation is theoretically established. The results of model predictions are 
assessed by comparison with available laboratory data in turbulent jet flows. 


Introduction 

Despite extensive recent contributions in direct and large 
eddy simulations of turbulent reacting flows, the application 
of such simulations is limited to “simple flows” (Givi, 1994). 
Based on this fact, it is now widely recognized that the “stat- 
istical” approach is still the most practical means in compu- 
tational turbulence, and future capabilities in predictions of 
engineering turbulent combustion systems depend on the 
extent of developments in statistical modeling. 

The literature on computational prediction of nonreactive 
turbulent transport is rich with schemes based on single-point 
statistical closures for moments up to the second order 
(Taulbee, 1989). Referred to as Reynolds stress models 
(RSM), these schemes are based on transport equations for 
the second-order velocity correlations and lead to determina- 
tion of “nonisotropic eddy-diffusivities.” This methodology is 
more advantageous than the more conventional models based 
on the Boussinesq approximations with isotropic eddy diffu- 
sivities (such as the k — e type of closures). However, the need 
for solving additional transport equations for the higher-order 
moments could potentially make RSM less attractive, espe- 
cially for practical applications. For example, it has been re- 
cently demonstrated (Hofler, 1993) that the computational 
requirement associated with RSM for predictions of three- 
dimensional (3-D) engineering flows is significantly higher 
than that required to implement the k - e model. The in- 
crease is naturally higher for second-order modeling of chem- 
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ically reacting flows due to the additional length and time 
scales that have to be considered (Toor, 1991; Jones, 1994; 
Libby and Williams, 1994). 

A modality to reduce the large number of equations associ- 
ated with RSM is to utilize “algebraic” closures. Such clo- 
sures are either derived directly from the RSM transport 
equations, or other types of representations (Speziale, 1991; 
Yoshizawa, 1988) that lead to anisotropic eddy diffusivities. 
One of the original contributions in the development of alge- 
braic Reynolds stress models (ARSM) is due to Rodi (1976). 
In this work, all the stresses are determined from a set of 
“implicit” equations that must be solved in an iterative man- 
ner. A somewhat similar method was applied to the heat-flux 
equation by Gibson and Launder (1976). Pope (1975) offers 
an improvement of the procedure by providing an “explicit” 
solution for the Reynolds stresses. This solution is generated 
by the use of the Cayley- Hamilton theorem, but is only ap- 
plicable for predictions of two-dimensional (2-D) mean flows. 
The extension of this formulation has been recently done by 
Taulbee (1992) and Gatski and Speziale (1993). In these ef- 
forts, the Cayley-Hamilton is used to generate explicit alge- 
braic Reynolds stress models that are valid in both 2-D and 
3-D flows. 

The objective of this work is to expand upon the formula- 
tion developed by Taulbee (1992) (also see Taulbee et al., 
1993) for predictions of turbulent flows involving scalar quan- 
tities (Brodkey, 1981). The specific objective is to provide ex- 
plicit algebraic relations for the turbulent flux of scalar vari- 
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ables. Both nonreacting and reacting flows are considered. In 
the latter, a second-order, irreversible chemical reaction of 
the type A + B -» P is considered in isothermal turbulent 
flows with initially segregated reactants (Brodkey, 1975; Toor, 
1975). The closure explicitly accounts for the influence of the 
Damkohler number and includes the mixing solution in the 
limit of zero Damkohler number. Similar to previous contri- 
butions, the starting equations are the currently available dif- 
ferential equations for the second-order moments. Accord- 
ingly, several previously suggested closures for the 
pressure -scalar gradients correlations are considered. The fi- 
nal results are compared with available experimental data in 
turbulent jet flows. 
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Theoretical Background 

With the convention that the angle brackets < ) represent 
the ensemble mean value of a transport variable and the 
prime denotes its fluctuations from the mean, the nondimen- 
sionalized averaged equations in space (x,, i * 1, 2, 3) and 
time (0 for incompressible, isothermal turbulent reacting 
flows are: 


with C €i - 1.44 and C' 2 = 1.92 — C <3 *. The parameter x rep- 
resents the correction for the dissipation equation in round 
jets and is given by x =((0/(f)) 3 {ftBS) 1 where {} denotes 
the trace, S represents the mean-flow strain-rate tensor, S tJ 
= [<?( w* )/<?*, + ^ 1/2, and ft denotes the mean flow 
rotation-rate tensor, -[aiu^/dxj f - a^u^/ax^/l. With 

C €j 0.89, the spreading rate of jet flows is correctly pre- 
dicted with the nonlinear stress-strain relation. 

Treatment of the scalar variable requires the solution of 
additional transport equations (Chakrabarti et al., 1995) for 
the reactants* covariance (Y a 'l£>, and dissipations 
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Here u i9 p, p, Y a9 Re , and Sc denote the ith component of 
the velocity vector, the pressure, fluid density, mass fraction 
of species a, the Reynolds number, and the Schmidt num- 
ber, respectively, while < <o a ) represents the rate of chemical 
reaction (( ui A ) = 

<<u a > = - Da((Y A ){Y B ) + (Y^)), (4) 


where Da is the Damkohler number. The algebraic formula- 
tion entails a two-equation scheme in which the Reynolds 
stresses and the scalar fluxes are expressed by nonlinear 
functions of the mean gradients and the time scales of the 
flow (Wang and Tarbell, 1993). The mechanical time scale is 
determined by the solution of transport equations for the tur- 
bulent kinetic energy ( k) = (m'w')/ 2 and for the turbulent 
dissipation 


<e) =«7 


\ * x , I 


For shear flows, these equations are (Pope, 1978): 


+ <<* tt yp+<^y a '>. (7) 

By neglecting the third-order mass-fraction correlations, the 
chemical-source terms in the expanded form read (no sum- 
mation on Greek indexes in all subsequent equations) 

< w a yp + < ^ y 0 ' > = - Da(« yjyj > + < y;y; >< y* > + « yj y; > 

+ Full resolution of the nonlinear interactions 

in the chemical scalar fields requires significant computa- 
tional effort in practical applications (Hill, 1976; Givi, 1989; 
Fox, 1996). The neglect of the higher-order scalar fluctua- 
tions for the configurations considered here is justified (Givi, 
1989), but cannot be recommended for general applications 
(Wang and Tarbell, 1993). In such applications, the single- 
point probability density function (PDF) or the joint PDF of 
the scalar variable provides the required information (Toor, 
1962; O’Brien, 1980; Dopazo, 1994; Fox, 19%). The inclusion 
of the PDF is not attempted here. 

There are several methods for evaluating the scalar covari- 
ance dissipation (Jones, 1994; Newman et al., 1981; Jones and 
Musonge, 1988; Borghi, 1990). By using the first-order term 
in the two-scale direct-interaction approximation, Yoshizawa 
(1988) develops a generic model for the scalar dissipation. An 
equivalent functional expression is obtained from Yoshizawa’s 
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results by making use of the time-scales ratios r Q = 2< k >< e a )/ 
«6><y„ » (hereinafter e a = c aa ) and replacing the diffusion 
effect term by the inherent gradient of the turbulent flux. 
The equivalent form of this equation including the effects of 
chemical reaction is (Adumitroaie, 1997): 
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in which the chemical-source term is of the form 

= - DaU ( < A „ > + < »< Y b > + « t Be > + < e „ s »< Y a >]. 

(9) 


To determine the magnitudes of the model constants the 
transport equation for r a as derived from Eqs. 5-8 can be 
used in the limiting case of mixing: 


<*> 1 dr a 
<«> r a dl 



+ 
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where the production terms are P = ~ {u f i u , j )d{u l )/dx } and 
P a = - (u'jYa)d(Y a )/dXj. In the experiments of Warhaft and 
Lumley (1978) on decaying heated-grid turbulent flows it has 
been observed that the magnitude of r a is in the range 0.6 ^ 
r a ^ 2.4. In the experiments of Beguier et al. (1978) on ther- 
mal turbulence in several thin shear flows it is indicated that 
r a = 2. Based on this information, using the procedure de- 
tailed by Jones and Musonge (1988) it is possible to estimate 
the magnitudes of the model constants: C y = r Q = 2.0, C = 
2.0, C, 4 « Q 2 - 1 - 0.92- C <3 *, C y2 = 0.5. ' 

To complete the closure formulation, all the third-order 
transport terms are described by the gradient diffusion hy- 
pothesis. Denoting by H any of the fluctuation products on 
which the second-order correlations rest, we have: 

<k) £?<H> 

< u ' h > = -c s — <«>;>— — , (id 

< € > ' dXj 


where C s is taken to be equal to 0.22 for all nongradient 
correlations (H = /c and 5 = 7 a ' 2 ), whereas for the dissipa- 
tions (H = e and H = e a ), C s = 0.18. Also, the molecular 
transport terms are neglected under the assumption of high 
Reynolds-Peclet numbers flow. 


Explicit Algebraic Models 

An improved explicit ARSM for 3-D flow has been derived 
by Taulbee (1992) from the modeled transport equation for 
the Reynolds stresses. This model is based on the general 
linear pressure -strain closure given by Launder et al. (1975). 
The improvement is due to an extended range of validity; the 
model is valid in both small and large mean strain fields and 
time scales of turbulence. The nonlinear stress -strain rela- 
tion for 3-D mean flows is of the form (Taulbee, 1992; Taulbee 
et al., 1993) a~a(S, H), where a is the anisotropic stress 
tensor, a i} ** [( w'w' )/(&) — 2 6, y/3]. The ARSM depends on 
key turbulence parameters such as the turbulence time scale 
t * (k)/(c); the product ion -to-dissip at ion ratio P/(e ), 
where P - -(k)a ij Sj i is the production of the turbulent ki- 
netic energy; the invariants of the strain rate and rotation 
rate tensors a 2 *(5^), and the model co- 

efficients of the pressure -strain correlation and the modeled 
dissipation equation. 

A similar line of reasoning is followed to obtain a 3-D alge- 
braic closure for the velocity-scalar correlations. The trans- 
port equations governing these correlations are transformed 
into algebraic expressions by making two assumptions: (1) ex- 
istence of a “near-asymptotic” state, and (2) the difference in 
the transport terms is negligible. The starting equations for 
the convective scalar fluxes are described by 


t a(u'K){ Uj ) «?(<»'»iX) + </> , y 0 ')/<p)a, J ) 
at axj dx t 

- Da((u'iY^XY p > + (u' i Y^XY a > + (u’Y^)) 



On the RHS of this equation, the following terms are iden- 
tified: turbulent transport, pressure -scalar gradient cor- 
relation, production by the mean velocity and the mean 
scalar gradients, chemical reaction effects, molecular trans- 
port (assumed negligible at high Peclet numbers), and viscous 
dissipation. Based on the Poisson equation satisfied by the 
pressure fluctuations one can arguably split the pressure- 
scalar gradient correlation into two parts corresponding to 
so-called rapid and slow terms (Lumley, 1978). The rapid 
term represents an inner product between the velocity gradi- 
ent tensor and a third-order tensor, the last one subject to 
symmetry, continuity, and normalization constraints. As sug- 
gested by Lumley (1978), since the slow pressure -scalar 
gradient term and the viscous dissipation term are functions 
only of turbulent quantities, they can be incorporated into a 
single closure. The ensemble of the entire pressure -gradient 
term and viscous-dissipation term enjoys a general relation 
encompassing some of the formulations proposed in prece- 
dent contributions. Consequently, this is written 
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The transport equation for the correlation coefficient ip ja (a 
# /3 ) is of the form: 
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The model coefficients in this equation are taken from Laun- 
der (1975): 

e la = 6.4, c, = 0.5, c, = 0; i = 2,6, (14) 

Jones and Musonge (1988): 

e 1 „-6/[l + 1.5(a y *a / *) V? ], c, = 1.09, c 2 = 0.51, 
c, = 0; i = 3,5; c 6 = 0.12, (15) 


Rogers et al. (1989): 



where Re, = 4(L) 2 /((e)j'), and Shih et al. (1990): 
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c, = 4/5, c 2 = - 1/5, c 3 = 1/10, 
c 4 = -3/10, c 5 = 1/5, c 6 = 0, (17) 


where H - 1.1 + 0.55(2i/i - 1) tanh [4(r„ - 1)]; i/» = l + 
F/18 exp ( — 1.11/ReY 1 )(12/ReY 1 + 80.1 log [1 + 15.6 (- // + 
1.15///)]), with F = 1 +27////8 + 9///4, a parameter involving 
the second invariant //= -1/2 and the third invariant 
III = - l/3a lj a jk a lc , of the Reynolds stress anisotropy tensor; 
Re t = 4<Jfe) 2 /(9(e)»') denotes the turbulence Reynolds num- 
ber: d jk = > - (u)u' k )(Y? »/ ((«;y:x«;y:>- 
2(Jt>(y^ 2 »; F D =9/2-27dj l /2 + 9d ?; // rf is the second in- 
variant of the tensor d ik \ djj = </ y/ d, y ; and dj, = d jl d, m d m) . 

To proceed, let us denote the mechanical -chemical corre- 
lation coefficient (normalized scalar flux) by: 
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where the notation D/D/ indicates the convective transport, 
and 7)“, 7} a , and denote turbulent transports of the scalar 
flux, the scalar variance, and the kinetic energy, respectively. 
Moreover P a ^ -{{k){Y^)) v2 if> ja d{Y a )/dx j is the produc- 
tion of scalar variance; S Q - {^ a Y^) is the chemical source 
term in the <K a ' 2 > equation; and the remaining quantities are 
the normalized production, pressure-gradient, and the chem- 
ical-source term: 
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Si. = -Da{v, a (Y p )+ <P ie (Y a ) + y lae {Yf) vl ). (22) 

Here y la „ = {u',Y^)A(k)(Y^)(Yp 1 )) V2 is the normalized 
covariance flux vector. 

The results of direct numerical simulations (DNS) of non- 
reacting passive scalar mixing in homogeneous turbulent shear 
flow (Rogers et al., 1989) suggest the existence of an asymp- 
totic state for the normalized correlation coefficient <p ia , but 
not for the scalar flux itself. This observation justifies the first 
assumption, at least for reacting flows near the frozen limit. 
The second approximation is yet to be substantiated and its 
assessment requires future DNS or laboratory experiments. 
Under these assumptions the term representing the convec- 
tive transport is set to zero and the difference in turbulence 
diffusion terms is discarded. This procedure leads to an alge- 
braic system of equations for the two unknown vectors ip ia 
and 


<Pa + D a^<P* + + C Q = 0 

<Pfi + + B n<P a + c 0 =0, 


1938 


August 1997 Vol. 43, No. 8 


AlChE Journal 



where the coefficients are 


the Cayley-Hamilton theorem and yields an expansion defin- 
ing a natural basis for this problem: 
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In the Appendix the inversion procedure via the Cayley- 
Hamilton theorem is outlined and the coefficients a n and a' n 
are listed. The final results provide an explicit solution for 
the scalar fluxes. In the limit Da -> « f the use of the mixing 
solution (Ofl^O) for the transport of a Shvab-Zel’dovich 
variable (Toor, 1962; Williams, 1985) is recommended. 
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Finally, the anisotropy of the turbulent diffusivity is ensured 
by the properties of the second-order tensor A: 

A ik = [(1 - c, - c 2 )S lk + (1 - Cj + c 2 )H lk - (c 3 + c A )a tj S jk 

- c 5 a k)Sji ~ (c 3 - c 4 )fl |; (30) 

This tensor turns out to be traceless (^,, = 0) as a conse- 
quence of incompressibility and of the particular values taken 
by the constants, c/s. Now, the solution of the system of Eqs. 
23 is conveniently represented in the form: 


j *>„ = -M-'[(« + D p A)C a - B a C p ] 
[ + D a A)C p - B e C a ), 


where M denotes the matrix [(1 - B a B p )6 + (D a + D^)A + 
D a Dp A 2 ]. The expressions for the turbulent fluxes of reacting 
scalars exhibit the influence of the Damkohler number Da. 
Also, the coupling between the reactants is reflected by the 
nonlinear dependence on the mean scalars and the presence 
of the covariance flux. 

To provide a computationally efficient algorithm, the ma- 
trix M is inverted analytically. This is achieved by the use of 


Illustrative Examples 

In this section, we present some sample results of numeri- 
cal calculations based on the models given earlier. There are 
two primary reasons for conducting these simulations: (1) 
model assessments via comparisons with laboratory data, (2) 
demonstration of the model capabilities in comparison to tra- 
ditional closures based on the linear gradient-diffusion ap- 
proximation. The flow configurations considered consist of 
turbulent-plane and round-jet flows for which laboratory data 
are available. The mean flow motion in these shear flows is 
assumed 2-D or axisymmetric. The space coordinates are 
identified by x -[*, y\ where x is the streamwise coordinate 
denoting the direction of the flow’s principal evolution, and y 
represents the cross-stream direction. The velocity field is 
identified by u«[u, u]. In nonreacting flow simulations the 
mass fraction of one conserved species, Y A , is considered. In 
the jet configurations, Y A * 1 is issued at the inlet into a sur- 
rounding of Y a = 0. For the reactive case, Y A - 1 is issued at 
the inlet into a surrounding of Y B = 1. These species are as- 
sumed thermodynamically identical, and there is no trace of 
one of these species at the feed of the other one; that is, 
complete initial segregation. Also, the heat generated by the 
reaction is assumed negligible. 

The transport equations governing the velocity and the 
scalar fields are of parabolic type with the thin-shear layer 
approximation. For 2-D mean flows, the ARSM (Taulbee, 
1992) is of the form:' 


a = -2C M r 


S+ b x gT(r 2 ^ -6 - 6 (2) j + b 2 gr(Sfl - fiS) 


(33) 


where 5 (2) = [6/ ; 2) ] * 1 for i — y * 1, 2 and 0 otherwise. The 
parameters and g are given by 

c W* 

1 “ -U>\gr) 2 <r 2 + 2(b 2 gr) 2 a> 2 




Ci + C -2 + (2-C )iy<€> + -^ 

a Dt 


(34) 


where C\, b lf and b 2 are constants from the pressure -strain 
correlation model (Cj *1.8, 6, » (5 - 9C 2 )/11, b 2 * (1 + 
7C 2 )/11, C 2 =*0.45). In the self-preserving regions of turbu- 
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lent shear flows, the convective term Da/Dt can be ne- 
glected. For 2-D mean flows, with zero rate of reaction the 
scalar-flux model is expressed as 


w>-- 


2rh a {k) d(Y a > 

1 + // G ^ dXj ' 


(35) 


which has the gradient form, but with an anisotropic diffusiv- 
ity. With the thin-shear layer approximation, 

II C * K C 2 c 4 fl 22^1 “ C l “ C 3 fl ll “ C 5 a 22^~ C 3 a 12^ 

( d(u )\ 2 

x K— ) <36 » 

and the nonzero components of the diffusivity tensor, A i; , 
are 


( d(u)\ 

2 ' 

A„-|l 2c 3 rh a a n ^ J 

(l-c 6 )a n + j 

d(u) 

2rh a (l c s )a 2l (1 c x c 3 a u 

dy 


2 ' 

A 22 “ 1 1 + 2 c 3 r/i a a l2 I 

(l-c 6 )a K + - 


4 22' 


d(u ) 

+ 2r/i a (l — c 6 )a 12 — — — (c 2 + c 4 a 22 ) (38) 


a 33 = (i + // c ) 


(l-c 6 )a 33 + T 


(39) 


d(u) 

A 12 “ I 1 — 2c 3 r/i a fl 12 | (1 “ c &) fl i2 


2 r/i a (l-c 6 )a„ + 




dy 


(1 - c, - c 3 a u - c 5 a 22 ) (40) 


A,, = l + 2 c 3 x/r a. 


d<«> 

*"^T 


(l“C 6 )fl 2 


+ 2 xA a |(l-c 6 )a n + - 


d<w> 




(c 2 + c 4 fl 22 ). (41) 


These anisotropic diffusivities are determined directly from 
the velocity gradient, the components of anisotropic Reynolds 
stress tensor, and the model coefficients. 

The numerical algorithm for the solution of the transport 
equations augmented by the algebraic closures is based on a 
first-order upwind differencing for the convection terms and 
a second-order central differencing scheme for all the other 
terms. Due to the anisotropic character of the algebraic clo- 
sures, it is possible to evaluate all the components of the 
Reynolds stress tensor and the scalar-flux vectors. In this 
evaluation, the terms appearing as model coefficients (e.g., 
P/(c) in Eq. 26) are treated in an iterative procedure. The 
implementation of the boundary conditions is similar to that 
in many previous simulations of parabolic shear flows (e.g., 
Taulbee, 1989). In the results presented below, the spatial 


coordinates are presented by 17 = yAx-jt 0 ) for hydrody- 
namic and 


y«^>-0.5<^>cJ 

for the scalar variables. x 0 denotes the virtual origin of the 
jets. In the nonreacting jets, the subscript CL denotes values 
at the center line (i.e., y * 0). In the reacting jets, the corre- 
sponding profile of Y A under no chemical reaction is em- 
ployed in the normalization. In all the figures below, the 
transverse variations of the statistical variables are presented. 

The experimental results pertaining to the velocity fields of 
planar jets as reported by Gutmark and Wygnanski (1976), 
Bradbury (1965), and Heskestad (1965) are compared with 
the model predictions in Figure 1. The agreement is reason- 
able for the mean streamwise velocity and also for the com- 
ponents of the Reynolds stress tensor. The predicted spread- 
ing rate ( dy <u)Au>CLWm0S /dx ) is 0.105, which is within the range 
suggested by experimental measurements. In Figure 2, the 
predicted results for the mean and the variance of the nonre- 
acting scalar are compared with the experimental data re- 
ported by Browne et al. (1984), Bashir and Uberoi (1975), 
Uberoi and Singh (1975), Jenkins and Goldschmidt (1973), 
and Antonia et al. (1983). Figure 2 indicates that the models 
based on the coefficients proposed by Launder (1975), Jones 
and Musonge (1988), and Shih et al. (1990) predict the mean 
scalar values in these experiments reasonably well. These 
models also yield good predictions of the experimental data 



Figure 1. Cross-stream variation of (u)/(u) C u <«' 2 )/ 
<u>ci, <v'*>/<u>ct, and <uV>/<o>c L for the 
planar jet. 
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of Uberoi and Singh (1975) for the scalar variance. All the 
other experimentally measured variance profiles are better 
predicted by the model with the coefficients of Rogers et al. 
(1989). 

The procedure by which the Reynolds stress tensor and the 
scalar-flux vector are determined by our explicit solution al- 
lows a direct comparison of the calculated fluxes with data. 
This comparison is made in Figure 2 and indicates that the 
models with coefficients of Launder (1975), Jones and Mu- 
songe (1988), and Shih et al. (1990) yield results in reason- 
able agreements with the experimental data of Jenkins (1976), 
but overpredict the experimental data of Antonia (1985) and 
Browne et al. (1984). These data are in better agreement with 
the predicted fluxes based on the model of Rogers et al. 
(1989). The lower spreading rates predicted by the model of 
Rogers et al. (1989) are primarily due to the relatively large 
values adopted by the parameter G la . In this model, the pro- 
posed form of C la and its correlation with the turbulence 
Reynolds number are determined with comparative assess- 
ments by DNS results of homogeneous turbulent shear flows. 
In the jet-flow experiments, as considered here, a direct ap- 
plication of the model yields relatively large values for C la , 
and thus small turbulent diffusivities. Consequently, the pre- 
dicted scalar spreading rate is lower than that measured ex- 
perimentally. Nevertheless, in the core region, the results 
predicted by this model are closer to the majority of available 
experimental data compared to predictions based on other 
models. 

With these results it is possible to perform an a posteriori 
appraisal of the closures based on conventional linear gradi- 
ent-diffusion hypotheses. For example, the parameters C M 
and Sc t as given by 



nar jet. 


<*V> 


Sy 


<*>%'> = - 


SCi 


v t 



dy ’ 


(42) 

(43) 


can be directly evaluated. The explicit algebraic relation for 
C M is given by Eq. 34; the relation for the turbulent Schmidt 
number is 


for all these parameters near the free stream. The amplitude 
of the parameters at the free streams can be controlled by 
modifications of the boundary conditions. An extract specifi- 
cation of these conditions requires inputs from laboratory 
measurements. 

Some of the influences of the chemical reaction on the 
scalar field in the turbulent plane jet are presented in Fig- 
ures 4 and 5. In the calculations pertaining to these figures, 
the model coefficients of Launder (1975) are employed. The 
influence of reaction in modifying the amplitudes of the 


Sc, 


1 -4[cf«i2 - (c 2 + c 4 a 22>(l - C, - 

d(u)\ 2 

C 3 a l\ C 5 a 22)]l^£i r 


fc a jjl+2c 3 T/. Bfll2 ^ J 

2' 

(l-cja^+j 

d(u ) \ X 

+ 2T/i„(l-c 6 )a 12 — (c 2 +c 4 a 22 )J 15 

l+(*2 2 


2g 


(44) 


Figure 3 shows the cross-stream variations of C u and of Sc t 
and r a based on the pressure -scalar gradient model in Shih 
et al. (1990). These results can be compared with - 0.09 
and Sc t = 0.7, typically employed in the linear gradient-diffu- 
sion approximations. Also, the ratio of the velocity to scalar 
time scales (r a ) indicates that an approximate constant value 
can be used at the central region of the layer. This is in ac- 
cord with the results of Beguier et al. (1978) and Tavoularis 
and Corrsin (1981). As expected, there are large variations 


scalars’ means (Figure 1), variances (not shown), and turbu- 
lent fluxes are captured by the model (Dutta and Tarbell, 
1989; Gao and O’Brien, 1991). In accord with the physics of 
turbulent flows with segregated reactants, the unmixedjiess is 
negative throughout the layer (Shenoy and Toor, 1989; 
Leonard and Hill, 1991; Wang and Tarbell, 1993). The same 
is true in the limit of no chemistry; in that case, the ampli- 
tude is slightly larger (Leonard and Hill, 1988; Frankel et al., 
1993, 1995). 
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Figure 3. Cross-stream variation of C M , Sc tJ and r* for 
the planar jet. 




Figure 4. Cross-stream variation of < Y a )/(Y a ) C l and 
(Y^Yg) for the nonreacting and reacting planar 
jets. 


The comparison between the full second-order model (Eq. 
12) and the algebraic closure in predicting the scalar fluxes is 
presented in Figure 5, and indicates that the transverse flux 
as predicted by the algebraic closure is in close agreement 
with that by the transport-equation model. However, there 
are differences between the two predictions of the stream- 
wise flux near the jet center line. The zero value of this 
flux in the algebraic model is due to the thin-shear-layer ap- 
proximation. The neglect of the axial diffusion in this 
approximation combined with the gradient diffusion nature 
of the algebraic model can only yield zero flux values at the 
axis of symmetry. While the thin-shear-layer approximation is 
also invoked in the transport equation model, the inclusion 
of the streamwise convective effects in the transport equa- 
tions can, and does, yield nonzero flux values. If the thin-layer 
assumption is relaxed in the algebraic model, the inclusion of 
axial scalar gradients would generate nonzero scalar-flux val- 
ues at the center line, thus reducing the disagreement. It must 
be noted that for this class of flows the cross-stream scalar 
flux is more dominant than the streamwise flux in influencing 
the mean scalar distribution and the production terms. Thus, 
the agreement observed in Figure 5 is encouraging in support 
of the algebraic approximation. Nevertheless, this is demon- 
strated here only for a “simple” flow configuration. The im- 



T1 



planar jets. 

It was computed via the algebraic scalar-flux model (ASFM) 
and the full second-order transport-equation model (SOM). 
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Figure 6. Cross-stream variation of <u)/<u>cl> <u' 2 V 
(w)cl» < v' 2 )/(u)cl , and < u'v')I(u)cl for the 
round jet. 


plementation of the model for complex flows would be very 
useful in further assessment of the algebraic approximation. 

The performance of the models for prediction of axisym- 
metric jet flows is assessed in Figures 6 and 7 where the ex- 
perimental data of Hussein et al. (1994), Abbiss et al. (1975), 
Wygnanski and Fiedler (1969), and Rodi (1975) are used for 
hydrodynamics variables, and those of Chevray and Tutu 
(1978), Becker et al. (1976), and Lockwood and Moneib (1980) 
for the scalar variables. The predicted hydrodynamic spread- 
ing rate with the axisymmetric correction is 0.094, and is in 
agreement with the experimental results of Hussein et al. 
(1994). Again, all the mean values are reasonably well pre- 
dicted. The same is true for the Reynolds stresses, except for 
the streamwise normal stress in the central region for which 
an improvement of about 30% can be obtained if all the com- 
ponents of the rate of deformation tensor are considered. 
Consistent with the planar-jet results, the model with coeffi- 
cients of Rogers et al. (1989) results in lower scalar diffusivi- 
ties. This model also yields lower values for the second-order 
moments in the core. The model predictions based on the 
coefficients of Launder (1975), Jones and Musonge (1988), 
and Shih et al. (1990) overpredict the experimentally mea- 
sured scalar’s covariance and turbulent fluxes. The predic- 
tions based on the model of Rogers et al. (1989) again yield 
better agreement for the scalar fluxes. It is important, how- 
ever, to indicate that the experiments of Chevray and Tutu 
(1978) are not conducted in the self-preserving regions of the 
jet. Therefore a definite assessment cannot be made without 
comparisons with further data. 



jet. 

From the preceding comparisons, it can be concluded that 
the algebraic model developed here provides an effective 
means of predicting the second-order moments in reacting 
turbulent flows. Because of their anisotropic feature, these 
algebraic schemes are more general than the conventional 
linear gradient-diffusion schemes (Toor, 1991). The explicit 
nature of the relations is particularly convenient for applica- 
tions in practical flows of the type considered by Hofler 
(1993). With the reasonable agreement of the model results 
with experimental data in simple configurations, the method- 
ology is recommended for predictions of more complex flows. 
Even so, the restrictions stemming from the assumptions in- 
volved in the development of the models have to be clearly 
underscored. The results shown here indicate the need for 
refinement of current pressure-gradient correlation closures. 
Moreover, the modeled transport equations for the passive 
scalar dissipations have known inconsistencies (Pope, 1983). 
These and the nature of the pressure-correlation models 
might raise realizability concerns, which can be considerably 
alleviated by implementing some of the techniques developed 
by Shih and Shabbir (1994). The present models are devised 
for high Reynolds-Peclet number flows; therefore some cor- 
rections might be required for modeling of the near-wall re- 
gions. In flows with important nonlocal effects such as the 
action of diffusion over long distances, rapidly varying flows, 
or other kind of flows far from equilibrium, the results by 
algebraic models are expected to have a larger departure from 
those by the full second-order moment formulation. In flows 
with very large strain fields there is a potential for singular 
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behavior of the scalar-flux models. This issue requires further 
investigations of complex flows. Further improvements are 
recommended by considering low Reynolds-Peclet number 
effects, the higher-order moments of the scalar-scalar fluctu- 
ations in reacting flows, and the effects of exothermicity in 
nonequilibrium chemically reacting systems. 
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Appendix 

The procedure leading to explicit solutions for the scalar- 
flux vector, as governed by Eq. 23 is described here. 

Consider an arbitrary 3-D second-order tensor Qm[Q.\ 
and the corresponding Kronecker tensor 8 = According 
to the Cayley- Hamilton theorem, this matrix satisfies its own 
characteristic polynomial: 

Q 3 -IqQ 2 + IIqQ-III q 8 = 0, (Al) 


(S + Gr '~M — [ G 2 + ( 2- WG 

JiI 6+G 

— /,+ c + //, +c )5]. (A4) 

It is easy to show that I t+C = I c +{«}, II I+G - 2 I c + II G + 
///*♦« ”Ie + He + Hie + ({5}/3). Therefore the nor- 
malized scalar flux vector takes the form: 


tp a =a 0 C + a,GC + a 2 G 2 C 
with the coefficients: 
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In a solenoidal velocity field the pressure -scalar gradient 
correlation includes a rapid term that satisfies the zero-diver- 
gence constraint. This further translates into {A) - 0; thus, 

1~{G 2 ) 

fl o * j j (A9) 

1 ~2 {C2}+ 3 {C3} 


1 , 1 . 
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where I Q = {Q) = Q u , II Q = 1/2 UQ) 2 - (Q 2 )] = 1/212,0 - 

QhQhI iHq = 1/6KC) 3 - m«2 2 ) + 2{<2 3 )] = l MQuQj’Qn 
-lQi,Q lk Q kj +2Q ij Q ik Q ki ] are the three tensorial invari- 
ants. Multiplying the characteristic polynomial with Q~ l and 
solving for the inverse, we obtain 

e -1 = 7^-(e 2 -/ G e + // 0 fi). <A2> 


1 

^2 1 } (All) 

1--{G 2 )+ -{G 3 } 


The reacting case is somewhat more complex. Nonetheless, 
by following the same procedure explicit solutions are 
obtained: 


+ * 0 C e + a,AC a + a\AC p + a 2 A 2 C a + a' 2 A 2 C B 
This relation can be used now to fmd explicit solutions to the p 

problem considered here. For example, for the case with Da (A12) 

= 0 (pure mixing), we can write 


Va - — (4 +G)~ l C, 


(A3) 


<Pp “ b o C a + b' 0 C p + + b\AC . + b 2 A 2 C a + b' 2 A 2 C B 


(A13) 


where G — D a A. Hence 


with the coefficients: 


AIChE Journal 


August 1997 Vol. 43, No. 8 


1945 



( [a*\ 

F. + D 2 ^ ) + B a B, | —jj—(DpF a - ED a )D, - £(£ + ^D a D, 


(l-B a B,)(F a F,-E 2 B a B,) 


■, (A14) 


F.(f, + D 2 ^) + D a ~(D a F, -ED,)- EB a B,^E + —^~D a D, 
' ( 1 - B a B,){^F a F, — E 2 B a B, ) 


, (A15) 


O a (F o F p -£ 2 B o B 0 ) 


/ U 2 > 1 B a B, \ U 3 

(A16) F„ = (1 — B a B l9 )l Z)„ — 


Bg D« F e ~ - E(D, - D a B a B,) 

D a F a F, - E 2 B a B, 


D„ F, - ED, B a B, 
° 2 ~ F a F,-E 2 B a B, 


D a F, ED, 
'F a F,-E 2 B a B,’ 


(A17) I {A 2 } 1 B„B,\ U 3 } 

Ji-u-*.v|n— - d c — dt )- d ‘— <a21> 

(A18) + (A22) 

The coefficients b, are obtained from the a- s through the 
(A 19) permutations a -» 6, p -* a, 0 -> /3, fl o “> ^o» *o ^o> fl i ~ * 

i>'„ fli -> i> lf fl 2 “* ^2* &2* 


with the shorthand notations: 
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I. Introduction 


There has been a resurgence of interest in effects of compressibility on turbulent flows related to 
the design of high-speed/high-altitude engines. Although experimental and numerical information 
is growing (for reviews see Refs. 1-3), rational theoretical and modeling efforts are in a preliminary 
stage of development. This is consistent with the fact that what is understood of the relative im- 
portance of many of the different physical effects of compressibility is very much in flux — changing 
as more numerical simulation information becomes available. While several issues regarding the 
effects of the compressibility of turbulent fluctuations have been recognized, progress in incorporat- 
ing the physics responsible for these new compressibility effects in single-point turbulence closures 
has been slow. One of the major impediments to progress has been the absence of a procedure 
that would allow for the inclusion of the effects of compressibility on the pressure-strain covariance 
appearing in the second-order moment equations. A method of including compressibility effects as 
they appear in the pressure-strain covariance and also the variable inertia effects are the subject of 
the present study. 

Several earlier studies have obtained closures for diverse effects of compressibility. Researchers 
have, in general, exploited a decomposition of the compressible field into solenoidal and dilatational 
parts. This has been done using a dimensional analysis in physical space 4 or in Fourier space, 5 
asymptotic analysis, 6 rapid-distortion theory, 7 and a singular perturbation method. 8 All such 
approaches have generated models for the scalar compressible terms, the pressure-dilatation and 
the dilatational dissipation. These scalar terms appear in the kinetic energy equation for high-speed 
flows. Such an approach to compressible turbulence modeling has been called, very sensibly, an 
“energetic” approach to the effects of compressibility. 9 The models resulting from the “energetic” 
approach to compressibility have been applied as compressibility corrections to the standard k — ( 
model, 10 and their generalizations, 5 or to standard incompressible second-order moment [Reynolds 
stress] closures. 11-13 Such a procedure implies, of course, a tacit assumption that compressible 
effects do not manifest themselves in either the pressure-strain or in the dissipation of enstrophy. 
Which is to say the effects of compressibility occur only in those terms explicitly linked to the 
dilatational field. Thus, energetic approaches to the problem of compressible turbulence, as pointed 
out in Simone et al., 9 are incomplete. 

At one time, the pressure-dilatation and the compressible dissipation, on which modeling effort 
has been expended, were believed to be the primary physical effects contributing to the reduced 
growth rate of the compressible mixing layer. Recent studies 14,15 have demonstrated that the 
dilatational effects on the mixing layer are, in fact, much smaller than once believed. In addition, 
more recent direct numerical simulations (DNS) suggest that the pressure-dilatation covariance is 
nominally more important than the compressible dissipation, contrary to early proposals. 4,6 The 
pressure-dilatation does not, however, account for the reduced growth of the mixing layer. It 
appears, as suggested in Ref. 16, that the phenomena responsible for the reduced growth rate of 
the turbulent shear flows is due to the reduction in the Reynolds shear stress anisotropy; this 
effect is thought to be due to the effects of compressibility on the pressure-strain covariance. This 
viewpoint is consistent with the earlier numerical studies. 17 The article 13 and the later 18 studied 
Reynolds stress closures in the context of the DNS of the homogeneous shear. The study showed 
that the inclusion of the [then] current compressible dissipation and pressure-dilatation models 
in the Reynolds stress turbulence closures led to improved predictions for the turbulent kinetic 
energy. However, there were no changes in the anisotropy consistent with that seen in DNS. The 
authors concluded that the [then] current models were deficient primarily in the modeling of the 
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pressure-strain covariance which controls the level of Reynolds stress anisotropy. 

It should be noted that the improved agreement for the time evolution of the kinetic energy, 13 ’ 18 
when using such models for the scalar compressible terms cannot be taken to indicate that such 
flows were rationally predicted. The models current at that time, designed with erroneous assump- 
tions regarding the importance of the dilatational effects, were providing dissipative behavior by a 
mechanism that did not reflect actual flow physics. This has since been substantiated numerically 
in the studies; 9,14 " 16 all of which indicate the lack of significance of both the pressure-dilatation 
and the compressible dissipation. The fact that the compressible dissipation and the pressure- 
dilatation are nominal effects is also consistent with the analytical development of Ristorcelli. 8 In 
Ref. 8 the pressure-dilatation is shown to vanish as turbulence approaches equilibrium; the sim- 
ulations mentioned are quasi-equilibrium flows for which the pressure dilatation is expected to 
be small. It should be mentioned that in the homogeneous shear simulations, arguably the most 
non-equilibrium of the benchmark flows, the pressure-dilatation is small but non-negligible; it is 
some 5-10% of the dissipation. Our position is the same as the position of Vreman et a/. 14 as 
pithily summarized in their conclusion. To paraphrase, turbulence models constructed using the 
dilatational dissipation or pressure-dilatation to explain the suppression of the turbulence [in the 
mixing layer] do not appear to be reflecting the correct physics. 

Thus, if it is assumed that the primary source of the reduced mixing rate in the mixing layer is 
due to the reduction in the shear stress anisotropy as indicated by DNS then a closure of the second- 
moment or Reynolds stress equations is in order: the pressure-strain covariance appearing in the 
second-moment equations is the only possible mechanism for such behavior. This is a substantially 
more difficult problem than that treated using the energetic approach; the quantities requiring 
closure are no longer scalars but second-order tensors. As one might expect, true compressible 
second-order modeling attempts are few. 7,19 

This article describes the development of a closure for the compressible aspects of pressure- 
strain covariance appearing in the Favre-Reynolds stress equations. [As is clear from the material, 
the phrase Reynolds stresses will be used to refer to the Favre-Reynolds stresses.] Closures for the 
unclosed terms involving the mean acceleration are also constructed. In as much as many engineer- 
ing calculations are done with lower order closures this article therefore also includes the additional 
development of the second-order moment closure into an algebraic Reynolds stress closure, used for 
flows near structural equilibrium, following established procedures. 20 " 28 

The algebraic Reynolds stress closure is noteworthy for the fact that it indicates that, even 
in the absence of mean deformation, the mean density gradient is a source of turbulence stresses 
in accelerating mean flows. As a consequence, for flows with large arbitrary mean density and 
pressure gradients an eddy viscosity representation for the Reynolds stresses is, from first principles, 
inappropriate. In the next section, Sec. II, the Favre averaged nondimensional form of the governing 
equations are given. Both first and second-order moment equations for a compressible medium, 
with no combustion, are given. Issues related to moment closures for compressible turbulence are 
also outlined. 

The development of closures for the effects of compressibility in the Reynolds stress equations 
is described in Sec. III. It is shown that, within the context of a pressure-strain closure linear in 
the Reynolds stresses, an expression for the pressure-dilatation covariance can be used to construct 
the off-diagonal components of the pressure-strain covariance. In this way the results of previous 
“energetic” approaches, 9 to the effects of compressibility can be built into the deviatoric portion of 
an expression for the pressure-strain. It is hoped that such a procedure would allow one to avoid the 
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development of a whole new theory and methodology for the compressible pressure-strain tensor. 

Closures for the effects of the mean acceleration, which involve the mass flux, are also developed. 
The new closures account for the influence of the turbulent Mach number, and the mean density and 
pressure gradients through a new quantity, the baroclinic dyad. The effects of the bulk dilatation are 
also included. Section 111 is concluded with a summary of models for the compressible dissipation; 
given the acknowledged lack of importance of the compressible dissipation in weakly compressible 
aerodynamic turbulence, as discussed above and in Sec. Ill, no development of models for the 
compressible dissipation is pursued. 

Starting with the closure for the second-order moment equations developed in Sec. Ill, Sec. IV 
develops an algebraic closure for the Reynolds stresses. The physics of compressibility, as captured 
by the full second-moment closure, are built into the simpler and more widely used two-equation k—e 
platform. The tensor polynomial representation techniques employed produce both two-dimensional 
and three-dimensional versions of an algebraic turbulence stress closure. Given the complexity of 
the three-dimensional algebraic closure and the current status of single-point turbulence models for 
three-dimensional flow only the two-dimensional model is developed into a working closure. 

Section V focuses on a numerical investigation of the closures. The numerical method used to 
simulate the free-shear flows of interest is sketched. The theory and results presented in earlier 
sections are implemented in simulations in the mixing layer. 

Simulations using second-order moment (SOM) closures as well as the algebraic models are con- 
ducted. The numerical experiments are constructed with the intention of investigating several very 
different issues of relevance to the prediction of compressible turbulent flows for engineering pur- 
poses. Foremost in importance is the success of an algebraic stress model to reflect the compressible 
physics of the full SOM closures. 

The pressure-strain methodology developed for the SOM equations in Sec. Ill are general and 
depend on a choice of closures for the pressure-dilatation. As consequence it follows that it is 
necessary to understand sensitivity of the formulation to different models for the pressure-dilatation. 
In this context two pressure-dilatation models are investigated. 

In addition to assessing the sensitivity of the pressure-strain model to different pressure-dilatation 
models and the suitability of the algebraic stress closure, we compare the computational results 
to what is seen in numerical and laboratory experiments. Of particular interest is the well known 
effect of compressibility on the reduction of the spread rate of the mixing layer. In as much as the 
reduction in the spread rate is due to changes in the principal axes of the Reynolds stress tensor 
the effects of compressibility on the anisotropy tensor are investigated. Given that the anisotropy 
of the Reynolds stresses is dependent on the pressure-strain, the effects of compressibility on the 
different components of the pressure-strain tensor are also investigated. 

II. Governing equations 

The problem formulation is now described. This includes a statement of the governing equations 
- the first and second-moment equations. Indications of the modeling issues to be addressed in 
subsequent sections are also given. The Favre averaging procedure is first described. 

The conservation equations for mass, momentum and energy in a Favre setting are now derived. 
The dependent variables used are the density p, the velocity u, and the total energy e t = h — 
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p/p + «,«,/ 2. The fluid is assumed to be Newtonian fluid satisfying the Stokes relation with zero 
bulk viscosity and a constant molecular Prandtl number. Although real gas effects are of interest 
for industrial applications there are a sufficient number of unresolved and more important issues 
associated with compressibility that justify limiting the study to ideal gases. Consequently, the 
pressure p is obtained from the equation of state p = pRT. 

The Favre averaging procedure is now described: denote by an over-bar the ensemble (or time) 
average and by the brackets the density weighted ensemble (or time) average: 

W = y (i) 

The ensemble average obeys the following decomposition rules: 

X = X + X\ X 1 = 0. (2) 

The Favre average obeys the following decomposition rules: 

X = (X) + X", (X") = 0, X 7i = X-(X). ( 3 ) 

The application of the above averaging procedure on the instantaneous transport equations, de- 
composing the variables Favre mean and fluctuating components, produces the Favre averaged 
equations. The equations of motion are first rewritten in nondimensional form (with respect to ref- 
erence values taken in the high speed stream: p <*,, Uqq, Too, /Xqq and the inlet value of the vorticity 
thickness i^). Using these reference quantities, we define the relevant nondimensional parameters: 
the Reynolds number Re = PooV-oo&u/ Poo, the Prandtl number Pr = c v p oo/A, the mean flow Mach 
number M = Uoo/VT^oo, 7 is c p /c v . For this study, Pr = 1. 


A. First moment equations 


The mass conservation equation reads: 

<TP . flp(ttj) n 

dt ^ d Xj 

and the conservation of momentum is: 


(4) 


&P{uj) . dp . dgji( u ) 

dt dxj dxj dx{ dij 

Note that the stress tensor notation is changed to <Tj,(u) = j£[S,j(ix) — 5Spp(u)<S,j] = 2 pS*j(u)/Re. 
For the strain rate Sjj(u) = + §|^) an d f° r all the other linear differential operators, this new 

notation is more suitable when investigating compressible flows. In these instances the two types of 
averages, having different properties with respect to the linear differential operators, are naturally 
encountered. Hereafter any tensor with a star superscript will indicate the deviator - the traceless 
portion of that tensor. 

Two supplementary hypotheses pertaining to the molecular transport of momentum are set 
forth: the viscosity fluctuations are unimportant and the mean viscosity (p) is described by the 
Maxwell-Rayleigh law, i.e., it varies with the mean temperature as {p)/p re j = ((T) /T re j) m , m = 
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0.76. The fluctuations of the viscosity and their correlation with other variables are neglected. 
Within the current notation the averaged stress is equal to <Tji((u)) + "). 

The Favre mean of the total energy of the fluid e t = T/[ 7(7 - 1)M 2 ] + njUj/2 obeys: 


&p{e t ) t &p{e t ){uj) 
dt dxj 


dxj 


dxj 


d 


( 6 ) 


As is the situation with the momentum equation, the molecular flux term in the energy equation 
introduces the other type of mean quantity present in the Favre formulation - the plain ensemble 
mean. Using the new notation, the heat flux 

m dT 

* (7 - \)RePrM 2 dxj 

has the resulting averaged expression: 

q j (T) = q J ((T)) + q j (T*). 

The remaining terms on the right-hand side of the energy equation can be expanded by using the 
equation of state in the average sense 


(7) 

Then the pressure work term is of the form 

PUj = p(uj) + 

the viscous work term reads 


u x a }i {u) = ( ui)aji((u )) + (u t )a J% {u") + u"a^((u)) + u" t <T_, t (u") 


and the turbulent total energy flux is 


*«?«?> = T(7 


The fluctuations of transported quantities in turbulent flows are sustained via interaction between 
the mean gradients and the turbulence. The fluctuations will be characterized by length and time 
scales of the order of those of the mean flow and which, provided the Reynolds number is large, 
will be many orders of magnitude larger than the fine scales at which the molecular diffusion is 
important. Thus the form of the moment equations carried in the closure development for the 
computation of free shear flows will not include viscous transport effects. 


6 



B. Second-order moment equations 


The average of the first moment of the Favre decomposed Navier-Stokes equations produces the 
following second moment equations: 

— ai — + — Bxl — = ~Wxl [p« u >fc) + p' u " s ik + 


+ * (w, + 1&) - <w, - “‘Z 


+U 


—dokiiiu)) , 


4 - u 

dxk ' dxk 


du" du" 


(8) 


It is necessary to provide a closure for two compressible quantities: the pressure-strain correlation 
and the mass flux/pressure gradient (the last three terms in the second line). The molecular 
diffusion terms are generally small in high Reynolds number flows and will be negligible for problems 
addressed here. It should be pointed out that, depending on how the derivation is done, the mass 
flux terms multiply the acceleration of the mean flow. For this reason the mass flux terms are often 
called acceleration terms: they appear to be important in accelerating mean flows. 

Various tensors are divided into their traces and their deviators. The production of the Reynolds 
stresses becomes, 


" - p » - 3 P4j - ~ p l ( “' + W ■ 

where P is the production of ( k ). The pressure-strain covariance is written as 




The viscous and pressure acceleration terms are written, respectively, as 

r kj( 

dx k 


V *- V 2 VS _ 777<^<« u » , 7T, d °kj(( u )) 2 -jy0*u((»» 

v \ 3 — 2 ““ U J ~r u i “ ~ u 


and 


= Mi, - -MSij = - 


dxk ' ~ x dik 3 _/ dxk 

-x&P , —jj dp 2-jjdp 


+ U‘ 


x dx, ' -’dxr^dxk"' 3 


(9) 


( 10 ) 


( 11 ) 


(12) 


Note that the mean pressure gradient appearing in M,j and the mean viscous stress appearing in 
V XJ can be replaced using the mean momentum equations. In which case they are written in terms 
of the mean flow Lagrangian acceleration. Both of these terms involve the mass flux; any closure 
for these two acceleration terms requires a closure for the mass flux. The dissipation is rewritten 
as 


du’! 


du " 


du" 

pZ 3 =p (f.j - h&a ) = °A u ")jrr + 


dx k 


'dx k 


dx k 


(13) 


In turbulent free shear flows the viscous diffusion part is overlooked owing to high Reynolds numbers 
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which are characteristic to these flows. The present analysis can accommodate the discarded viscous 
terms when necessary; the tensor will be carried for generality. Furthermore, the turbulent 
transport terms are considered together: 


T,jk = - [p(u"u"u'k) + j/u'jSik + p'u , -5j k - u'!o k j{u") - . 


The Reynolds stress or second-order moment equations are then written 

dp(u"u") + _ _d 


dt 


^ + ^ + n « + "*■ + *5 - ? * 

+| {p + Jd+M + V-pe]^ 


(14) 


(15) 


The equation for the anisotropy The single-point anisotropy tensor, a, is defined 

a „ = <»"<>/« - | Si,. (16) 

The anisotropy has zero trace. Note that the Reynolds stress has been normalized by { k ) and not 
2 (k); in which case this version of the anisotropy is related to an older definition by = 2b jj . The 
turbulence time scale is r ~ {k)/c s where e s is the solenoidal dissipation. The second invariants of 
the mean strain and rotation are denoted by a 2 = S£((u))5£((u)) and u 2 = Qij((u))Qij((u)). 

In order to derive the algebraic stress closure in Sec. IV the second-order moment transport 
equations are replaced with the equation for (k) and an equation for the anisotropy tensor. For 
the algebraic stress closure the modeling procedure of Taulbee, 21 discussed further in Sec. IV, 
is employed. Taulbee’s procedure, in order to be consistent with observed asymptotic behavior, 
chooses as dependent variable the ratio a XJ f(ra). The relevant form of the anisotropy equation is 


rap 


D a l} /(ra) 


Dt 



1 

<*> 


'dT iik 

(«> dT k 

Oij r dT k art 1 

dx k 

(k) dx k 

(k) 1 dx k dx k \ 

1 

+ W) 

[pr } + + m*j + v* - p f* j - 

-C„)-y 
P *3 

T Do n M + V + pfd-pj c 
+ o Dt + pe t 


(17) 


where D/Dt indicates the mean Lagrangian derivative, D/Dt = ft + ( u i)fr 3 - 


The kinetic energy equation 

dm wxtij) 

dt dxj 


The equation for the turbulent kinetic energy ( k ) = 

^ [?«*> 

-r, dp —doiiUu)) ~ ~M 

+ ^-"?4 +u; 


(«)/2 is 


(18) 


The additional terms reflecting the compressible nature of the turbulence are on the last line; 
they are, respectively, the pressure-dilatation covariance, the mass flux - mean acceleration, and 
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the compressible dissipation. Note that neither the pressure-strain nor the mass flux/pressure 
acceleration terms, M ,j, appear in the { k ) equation. The effects of the pressure-strain and the 
Mij , appears only in the Reynolds stress equations. Any classical two-equation turbulence model 
cannot, as a consequence, account for the physics associated with these two unknown covariances 
that lead to changes in the Reynolds stress structure. Early approaches accounting for the effects 
of compressibility have focussed only on the effects of compressibility as they occur in the energy 
equation. This so-called “energetic” approach, to use the phrase of Simone et a/., 9 for the effects of 
compressibility misses the changes in the Reynolds stresses due to compressibility and inertia. It is 
for this reason that, in Sec. Ill, the second-order equations are closed. An algebraic Reynolds stress 
closure is then derived in Sec. IV. In this way the structural effects of compressibility, for a certain 
class of flows, can be accounted for within the context of a two-equation single-point closure. 


The dissipation equation The dissipation appearing in the kinetic energy equation is typically 
written, for locally homogeneous flows, 4,6,29 


CTj, ( U ")fe = p € = p ( € ’ + €c ) = + \pSlk( u ") 


(19) 


with Slij( u) = as the rotation rate. It is customary to compute the solenoidal dissipa- 

tion from the incompressible k — t model extended to variable density flows: 


dp is dp d_ „ _ p (k, 

dt dxj dxj p ' u i €3 ' Redxj 




0 <*> 

dij 



( 20 ) 


where C { , = 1.44 and C <2 = 1.92. 


In the present closure formulation turbulent transport of the second order quantities are all 
closed using a gradient diffusion model. Thus the turbulent transport of the product of fluctuating 
quantities, S, is written 


«' s > = - c -f ■ 


( 21 ) 


where C s is taken to be equal to 0.22 for all non-gradient correlations (e.g., S = (k) or S = (u'£u f {)) y 
whereas for the dissipation (S = ? 5 ), C s = 0.18. 


III. Compressible closures for the second-order moment equations 

This section describes the development of closures for the unknown terms in the second-order 
equations. The pressure-strain covariance appearing in the Reynolds stress equations is closed using 
a linear tensor polynomial representation in the Reynolds stress following well established proce- 
dures similar to Launder et ai 30 In recognition of the minor role viscosity plays in high Reynolds 
number free shear flows, only pressure-strain and pressure- acceleration closures are addressed. The 
pressure-acceleration terms are closed using a leading order [isotropic] gradient transport expression 
for the mass flux. The following section, Sec. IV, uses the results of the present section to obtain 
an algebraic closure for the Reynolds stresses. 
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A. A closure for the pressure-strain covariance 

In incompressible turbulence the closure of the pressure-strain covariance by a tensor polynomial 
linear in the Reynolds stresses is used with notable success for simple two-dimensional mean flows. 
The procedure is standard; a popular early reference is Launder et a/. 30 An updated version of 
the linear 30 rapid pressure-strain modeling is given in Speziale et al 31 As the ultimate goal is 
to devise an algebraic Reynolds stress model only pressure-strain models linear in the Reynolds 
stresses will be considered. A discussion of nonlinear rapid- pressure-strain model and the need for 
realizability can be found in Ref. 32. Additional discussion of the physical assumptions underlying 
such methodologies can be found in Refs. 33-35. 

The compressible correction to the pressure-strain covariance representation is obtained using es- 
tablished linear procedures. As in the incompressible situation, 30,34,32 the pressure strain-covariance 
closure can be written as: 

“ P = Aij + 2 P [Zpiqj + Zp]qi][Spq({ u )) + ^p?(( W ))]* (22) 

Such an expression is possible if the supplementary compressible terms appearing in the Poisson 
equation for the pressure (see Appendix) are, in the weakly compressible limit, of higher order. 
This will be the case if the evanescent wave portion of the initial value problem is not important 
as is the case for aerodynamic applications. 36 The tensors in the above decomposition are modeled 
as: 


Ai 

— C\p 6a , j + A P p g 

(23) 

T 

J'pjqt 

(k) 

— fiqi^pj d“ ^2^pq^ij “I" ^qj^pi) d" Pl^pj^qi 

“h 02{&pqQ , ij d" &p\Q>qj “f* &ij &pq d - ^jq^pi) d“ i^^qi^pj 

(24) 


Here Zpj 9 », following precedent, is linear in the Reynolds stresses satisfying necessary symmetry 
requirements. There are five unknowns. To determine the coefficients in T p j q i additional constraints 
are required. As is the usual procedure, the normalization constraint, 30,32 requires 

Zppqi = [(3 + 2 a 2 )S qt + (3ft + 4f3 2 )a qi ](k) = «u"> (25) 

which can be satisfied if 


3 a x + 2a 2 = 2/3 
3 /?i + 4/3 2 = 1. 

There are now three unknowns and additional information is required to obtain them. 

In incompressible turbulence modeling the trace of the pressure-strain is zero, T p k q k = 0, and 
this provides the additional information to determine the unknowns. In compressible turbulence 
the trace of the pressure-strain is the pressure-dilatation and the so-called continuity constraint 
becomes 


A k k + Ap l pkqk [S pq + ft P ,] = 2 ptd. (26) 

As the right hand side, pd, is known from earlier energetic approaches to the compressible turbulence 
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modeling problem, the continuity constraint, Eq. (26), becomes a constraint that determines the 
coefficients in the pressure-strain closure. Note that as pd vanishes with turbulent Mach number 
the incompressible limit, Akk + 4 ~p X p kqk[S pq + = 0, is recovered for vanishing compressibility. 

As a consequence of the continuity constraint a certain combination of coefficients appears in 
the final model. These are readily defined by a portion of Eq. (26) such that 


Zpiqi/ (k) — (c*l + 4a2 )S pq + (/ 3\ + 5/? 2 + 03 ) a pq “ d\S pq + ^2^* (27) 

The final model for the compressible portions of the pressure-strain is then written in terms of d\ 
and c ?2 • The values of d\ and d 2 , as will be indicated shortly, are then determined by the expression 
for the pressure-dilatation. The values of the a, can be related to the d t : 


Qi = -dx/ 5 + 4/15, 
0i = (15 + 6C 2 )/33, 
/? 3 = d 2 + C2/2. 


a 2 = 3di/10 — 1/15 
02 = -(2 + 3C 2 )/22 


and C 3 = (5 - 9C 2 )/11, C 4 = (1 + 7C 2 )/11. 

Application of the normalization and continuity constraints then allows the linear pressure- 
strain model to be written as 


n* - p f *. = -ccp lay + m [(| + f^i) 5* ({«))+ 

[1 - C 3 + 2 d 2 ] [^.((u)) + S* p ((u))a p j - ^ p *,({ u »a p ,^] - 
[I-C 4 - 2d 2 ][a tp n pj ((u)) - n ip ((u))a pj ] + ^d 2 5 pp ««))a ti j 


(28) 


which follows from Eq. (22) after subtracting the pressure-dilatation from the left side and the 
pressure-dilatation model from the right side of Eq. (22). All terms involving the d, represent 
corrections due to the compressibility of the fluctuations; the di vanish as the turbulent Mach 
number vanishes. The terms involving the C x come from the incompressible pressure-strain model. 
In which case the choice of the C, allow one’s favorite incompressible pressure-strain model to be 
used. The above expression for the pressure-strain, reflecting the choice C 2 = 0.45, is used in all 
the calculations presented in this paper. In more complex flows realizability corrections may need 
to be incorporated, in which case C 2 would be a variable. A realizable form of the model is given 
in an Appendix; more details can be found in Ref. 37. 


1. Commentary 

It has been seen that the knowledge of one invariant of the pressure-strain, the pressure- 
dilatation, and the assumption of a form linear in the Reynolds stress allows one to obtain a 
model for the deviatoric components of the pressure-strain tensor. In this way an independent new 
theory for the compressible pressure-strain can be avoided by using the results of developments 
for the scalar pressure-dilatation. The results of previous so-called energetic approaches, Simone 
et a/., 9 to the effects of compressibility is built into the deviatoric portion of an expression for the 
pressure-strain. 
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This point merits consideration from another point of view. Consider, for the moment, the 
following partition of the pressure-strain: 

2 

psij = psr + -pdSij. (29) 

For this subsection, the primes on p have been dropped. Here by s*j we now mean the deviatoric 
portion of the fluctuating strain, s*j = 0, which contains solenoidal and compressible contributions 
s*j — s*f -hsj'j’ 7 , where, of course, s *j — sjj . This is done for ease of presentation; in the nomenclature 
of Sec. II these quantities would, of course, be represented by a precision unnecessary for 

the present discussion. Thus one can write 


P$ij = ps io 


+ p s ij c + \p d6 ' ij- 


(30) 


The term ps*- 1 is closed using standard incompressible pressure-strain closures. The pressure- 
dilatation is cl osed using models already in the literature, see Sec. III.B below. In Sec. III. A an 
expression for ps has been obtained. 

As has already been discussed in energetic approaches, the pressure-dilatation cannot account 
for the suppression of the turbulence by the reduction in the shear stress as is seen in compressible 
flows. A straightforward extension of the energetic approach to the second-order closure level 
implies 


PSij = P$ij 1 


2— i r 

+ 3 pdSij . 


( 31 ) 


Such an expression, as is consistent with Refs. 13,18, is likewise unsuccessful in reducing the turbu- 
lence shear stress. As will be discussed further in Sec. V the pressure-dilatation is a small quantity 
and has a nominal effect, for a a* 1, on the turbulence stresses and energy when it is included in 
the spherical portion of the pressure-strain. Our procedure, using the pressure-dilatation, produces 
an expression for the compressible contribution to the deviatoric portions of the pressure strain, 

Q 

ps *■ . In fact, as was born out by computational experiments, had the expression 


PSij = PS* XJ + PS 


u 


(32) 


[without pd on the diagonal] been used in the numerical investigations reported in Sec. V, our 
results would not have changed much, for a < 1. 


2. Pressure-dilatation models 

To obtain the final form of the pressure-strain an expression for pd is required. There are some 
choices. Some proposals for pd require the solution of transport equations for the density variance 5 
or pressure variance. 7 There are two pressure-dilatation models 8,38 that do not require separate 
equations. The model of Sarkar 38 is: 

In pp = 73 = -***? l|(^ - + ^ U <**> 
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where \\ — 0.15, X 2 = 0.2, and xz = 0 (still to be determined by the author). Note that the model 
as we have grouped the terms makes it appear as if it is singular in M t ; it is most definitely not 
singular. The model of Ristorcelli 8 is 


p'd = - p 1 + T k - -M?y(7 ~ 1 )(Pt + P ? + T t )] 

— / / \ n*2 /D(3<T 2 + 5u> 2 ) 

-p(k)Mf X ^ L 


where 


X = 


2 I vd 


1 + 2I pd M? + llpdMfrft -1)’ 

I P d=\l[ + r pd [ 3 tr 2 + 5 u; 


X = 




1 + 2/^2 + §7^7(7 -1)’ 

^ = s(I) 3 ° 2,r - 


(34) 


(35) 


Here a is the proportionality constant in the Kolmogorov scaling; we return to this in more detail 
below. In the above expression, Tt is the transport of the mean temperature including effects such 
as heat flux and the turbulent or pressure transport. The production term Pp, contains the mean 
pressure dilatation and the mean dissipation with positive sign (heat release term). These two 
models, as well as a third (the Aupoix model), have been discussed in Ref. 39. 

The d{ are determined by the pressure-dilatation model. For the Sarkar model, one requires 
that 


8x3 M? xiMt 

“1 = — 7 — > “2 = —77— 


' PP 


= XiMfp e 3 . 


For the Ristorcelli model one requires that 

xM t 2 


d 1 = d 2 = 


XM? 


(36) 

(37) 

(38) 


= xM 2 [p I + |m? 7 ( 7 - l)(Pr + pe)]~ p<fc)Af t V £>(3<T ^' 5h,a) - (39) 

Transport terms have been neglected; this is necessary to obtain an algebraic closure [which requires 
homogeneity]. The compressibility of the turbulence will only be important where the turbulence 
energy is large [and thus also M<], which is typically in regions of large production and where 
transport is not as important. Which is to say that transport is only of importance in the periph- 
eral regions of simple flows - regions where the turbulence intensity is low, production low, and 
[therefore] the fluctuations are essentially incompressible. 
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B. The mass flux and pressure-acceleration closures 


The models for mass fluxes have not undergone much development and their importance in a 
general flow is not fully understood. Judging from the equations one can infer that the mass flux 
will be important in accelerating flows with large density gradients. There has been some research 
on the mass flux: Taulbee and VanOsdol 5 solved its transport equations and Ristorcelli 40 showed 
how an algebraic model for the mass flux can be obtained from its transport equation: 


P u" = r u 


r d(ui) 2 d(ui) diuk) 

+ v 2 tI- X x ' ' k/ 


dij 


dxk dij 




(40) 


where r u = M*r/[ 1 + ~^(P / (p c) - 1)]. Here v o, u\ and v 2 are the coefficients retrieved from 
the inversion of the matrix Gij = ^0 = -(1 + Ig + //g)*^ ^ = (1 + I G )u 2 , 

V 2 = (1 + /g + //g + 77/g) _ 1 , the Roman numbers representing the invariants of G. Note that the 
leading order term is a gradient transport model. For an isotropic turbulence, an eddy viscosity 
formulation is possible. For the present set of simple benchmark shear flows, one does not expect 
the mass flux terms to be very important - the mean flow does not accelerate much. We shall use, 
for the sake of simplicity, the eddy viscosity form of the mass flux - the lowest order contribution in 
the polynomial given above. This is consistent with the gradient transport model used in Ref. 11. 

The baroclinic dyad, a tensor product formed from the mean density and pressure gradients, is 
defined as 


TZ{j — 


dp dp 
dx{ dxj 


(41) 


and the mass flux/acceleration terms in the second-moment equations can then be written, using 
the leading order portion of Eq. (40) - equivalent to replacing the Reynolds stress by two thirds of 
turbulent kinetic energy times the Kronecker tensor, as 


^ 1 ;— _ g (^ 1 ; + 7^71 g7£pp$,j) — g r uVo(k)Rij 

where R * = TZ tJ + 1 l Jt - §7 Z PP 6 XJ . 


(42) 


C. The compressible dissipation 

There are several models available for the compressible dissipation. 4-6,8 Many of these models 
reflect certain assumptions regarding the importance of the compressible dissipation observed in 
early DNS of compressible turbulence. The compressible dissipation has since been found to be less 
important than originally believed. In fact, the low turbulent Mach number asymptotics 8 indicate 
it varies inversely with the Reynolds number and as a consequence is negligible in engineering flows, 
though perhaps not in low Reynolds number DNS. Compressible dissipation models are nonetheless 
included for completeness. 

As the Taulbee and VanOsdol 5 compressible dissipation model requires additional transport 
equations, in the spirit of computational simplicity, that model will not be used. The closure 
proposed by Sarkar et al 6 is based on ideas from linear acoustics and appears related to the initial 
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value problem. 36 It can be written as 


? c = 3a s e„M? (43) 

with a a — 1.0 from DNS of decaying compressible turbulence; unfortunately this arrangement was 
deemed to lack universality. 19 Zeman 4 provides a model on the grounds that eddy shocklets occur 
in high speed flow and relating this assumption to the dilatational dissipation: 

£ C = 0.75(1 - exp {-[(|(1 + 7 )M t - 0.1)/0.6] 2 })e, (44) 

It has been shown by Blaisdell et al. 41 that Zeman ’s model gives incorrect scaling between e s and 
£ c . Besides, the exponential dependence on M 2 delivers a steeper growth rate reduction compared 
to other models. 

As has been mentioned in more recent DNS studies, 9,14,15,42 have demonstrated that dilatational 
covariance closures with a M f scaling predict effects of compressibility when they are, in fact, very 
small. In an asymptotic analysis, Ristorcelli 8 has found that the compressible dissipation has an 
M 4 dependence and is inversely proportional to the turbulent Reynolds number: 

?c = { 3o| ^2 + MV3] + (5) 7-2 [ 3ff2 + 5<i ' 2 ] • 

( 45 ) 

The parameters are Jf = 0.3, 7| = 13.768, 7J = 2.623, 7f = 1.392, 7£ = 3 and a r = 0.4-4 is the 
Kolmogorov scaling coefficient. Also, A7j denotes the turbulent Mach number, that is M 2 = §(fc)/c 2 
and R t the turbulent Reynolds number R t = p4(k) 2 /(9e s fi)Re. The local speed of sound is given 
by c 2 = T/M 2 . For high Reynolds number flows, in the absence of wall effects, viscous diffusion is 
negligible. 

As this article focuses on a closure for the Reynolds stresses, the turbulent temperature flux 
and its transport equation modeling is not presented here. The reader is referred to the thesis of 
Adumitroaie 37 for a treatment of the thermodynamic modeling issues. It should be made clear that 
the findings, for the classes of flows with which we are concerned with in this article, are insensitive 
to the turbulent temperature flux modeling. For simple unidirectional shear flows without strong 
heat transfer effects the difference between a second-order closure and an eddy viscosity closure is 
negligible. This cannot be expected to be true for more complex flows. 

At this point the Reynolds stress equations have been closed and it is possible to compute 
the flow using a second-order closure. Such computations are the subject of Sec. V. There are, 
however, a wide class of engineering flows that can be computed using simple algebraic Reynolds 
stress models. An algebraic closure for the Reynolds stresses is now developed. 

IV. A compressible algebraic Reynolds stress model 

A quasi-explicit algebraic model for the Reynolds stresses is now derived. An algebraic Reynolds 
stress model comes from the fixed point solution of the evolution equations for the anisotropy tensor. 
These equations can be thought of as describing a turbulence in a state of structural equilibrium: 
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the fixed point solution corresponds to an exact solution of the Reynolds stress equations. 

Several permutations of quasi-explicit algebraic Reynolds stress expressions exist. 2C> ‘ 28 The 
qualifier “quasi-explicit” is used to indicate that, as the fixed point equations are nonlinear, the 
solution is given implicitly. A notable exception is the recent explicit algebraic model of Girimaji 28 
who has found the exact nonlinear solution to the fixed point equations. The inception of our 
work predates 28 and our procedure follows precedents set by, 20 * 21 * 26 The polynomial representation 
methods will be used to obtain two-dimensional and three-dimensional versions of the algebraic 
turbulent stress models. 

The method of algebraic stress modeling introduced by Taulbee 21 is used. This involves the 
ansatz 


— [— 1 = 0 
Dt tct 


(46) 


which allows relaxation effects to be built into the Reynolds stress model. Equation (46) allows 
a relaxation of the anisotropy to its equilibrium value at the same rate that the relative strain 
reaches its equilibrium value. 21 The major improvement due to the formulation in Eq. (46) comes 
at small applied strains. It can be shown that Eq. (46) produces a stress model consistent with the 
weak strain expansion of the Reynolds stress equations; the usual algebraic stress modeling ansatz y 
jy t a,y = 0 does not. 

The combination of transport terms in Eq. (17) is set to zero following established algebraic 
stress modeling procedures: 


1 

'dT ijk 

(«) dT k 

a ij 

\dT k 

T m i 

(k) 

dx k 

(k) dx k 

(k) 

dx k 

8x*J 


(47) 


Applying these approximations results in a quasi-linear tensor expression for the anisotropy. 

If == | — 62 = =p\t u v 0, 63 = C3 — 2^2, 64 = C\ + 2d 2 , and 


C X J + C t2 - 2 + (2 - C £l )^- + + y (1 - 2 d 2 )S pp ((u)) + 2 — — pIc 

£5 p € s o Ut 6 p e 3 

the algebraic fixed point form of the Reynolds stress anisotropy equation is written 


i -1 


a = —gr 


b\S* + 62 R* + b 3 (aS* + S*a - |{aS*}<s) - 6 4 (afi - fia)] 


(48) 


(49) 


where the curly braces signify the trace. From this expression it is seen that the anisotropy tensor 


dij = a tj ( S*,fi,R*) 


is dependent on three second order tensors, two symmetric and one skew-symmetric. 


(50) 


A. A three-dimensional algebraic Reynolds stress model 

Standard representation theory methods can be applied to obtain the solution a,y = a^( S*, fi, R*) 
to Eq. (49). In contrast to the incompressible case the solution of Eq. (49) is now much more dif- 
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ficult inasmuch as the procedure now involves an additional tensor. Following standard methods 
the solution can be expressed as a finite 3-D tensor polynomial, 

a = ^C X T X (51) 

A 

comprised of a linear combination of all the independent tensor products (generators) formed from 
the three primary tensors. The coefficients in the polynomial are functions of the independent 
invariants of the tensors. For this problem the dimension of the [minimal] tensor base is A = 41 
(cf. Spencer 43 ). This is very large and is unlikely to be used in practice. 

The complexity presented by such a large tensor basis can be side stepped by simplifying 
approximation regarding the 63 term. It has been argued 21,27 that for the range of values used for 
the constant C2 the inequality C3 <C C 4 holds and thus the term multiplied by C3 will only have a 
small effect on the solution. This approximation decouples the contributions of S* and R* to a. As 
the equation is linear the solution is determined using the superposition principle. This one allows 
to split the problem into two equations of lower tensor base dimension: 

a = a s + a R (52) 

where a s stands for the solution dependent on S*, 

a s = -gr [b x S* - 64 (a 5 ft - fta s )] (53) 

and a R denoting the solution dependent on R* 

a fi = -gr [& 2 R* - 64 ( a R fl - fta R )j . (54) 

The decomposition of a into portions dependent on S* and R* is unique. 

Applying the results from Ref. 21 the strain dependent portion of the solution for the anisotropy 
tensor can be written 

a s = — 2CVrS* - 4a 2 r 2 (S*ft - ftS*) - 8a 3 r 3 (ft 2 S* + S*ft 2 - -{S*ft 2 }<5) 

3 

-16a 4 r 4 (ftS*ft 2 - ft 2 S*ft ) - 32a 5 r 5 {S*n 2 }(n 2 - i{ft 2 }<$) (55) 

o 

where = 615(1 + \hlu) 2 )h u a 2 = ^b i b 4 g‘ 2 h 2 , a 3 = <*4 = -f&i b%g 4 hi, a 5 = lhb^g s h u 

ho = h 4 gr, hi — /i 2 [l + 2/ijju; 2 ] -1 and h 2 = [2 + /iqU> 2 ] -1 . 

Using a similar procedure the portion of the Reynolds stress anisotropy dependent on the 
baroclinic dyad is written 

a R = — 2C M rR* - 4a 2 r 2 (R*D - ftR*) - 8a 3 r 3 (ft 2 R* + R*fi 2 - ~{R*ft 2 }<5) 

u 

-16« 4 r 4 (fiR*n 2 - ft 2 R*ft) - 32a 5 'r 5 {R*ft 2 }(ft 2 - -{&}&) (56) 

o 

in which the coefficients have the same form as those in a s with the exception that 61 is replaced 
by b 2 . To obtain the full anisotropy tensor the two complex expressions Eq. (55) and Eq. (56) need 
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to be added. In the light of the complexity of the three-dimensional formalism a two-dimensional 
formalism is developed. 


B. A two-dimensional algebraic Reynolds stress model 

A simpler and tractable two-dimensional treatment is possible. In many engineering flows the 
mean flow and the statistics of the turbulence are two-dimensional. The two-dimensional problem is 
less complicated as the number of tensor products necessary to express the solution is substantially 
reduced. The symbols S,fi,R now denote two-dimensional tensors. 

It is necessary to recast the equation for the anisotropy in terms of the traceless 2-D tensors: 
±Uj(( u )) = Sij((u))- ls pp ((u))slj\ igj(p,p) = Rij (p, v)-\ Rpp (p, p) 8$ . Here, the two-dimensional 

Kronecker symbol is = [<S f ^] = 1 for i = j = 1,2 and 0 otherwise. The pressure-strain model 
is then written: 

n* - p e* = -C{p ea ij + p(k) [[| + !*] S* «u»+ 

[1 - C 3 + 2 d 2 \ [a.pS^ti)) + S^ p ((u))a p j - |5;,«t t ))a p ^ li ] - 
[1 -C 4 - 2d 2 ][a ip n pj ((u)) - D,p«u»a w ]+ 



The pressure acceleration is written 


- — + ~T u V 0 (k)R pp 


The fact that both 2D and 3D expressions of the pressure-strain correlation model must give 
the same result when applied to two-dimensional mean flows will be used for our simplifications. 
Recasting the model in 2D is done to take advantage of the simplifications that result from the 2D 
structure. 

Inserting the closures for the pressure-strain and mean acceleration terms, Eqs. (57), (58) into 
the Reynolds stress equation, using the same ansatz regarding transport at a fixed point, the 2D 
analog of Eq. (49) is obtained: 


= -QT 6iS* + 6 2 R* + b 3 (as* + S*a - |{aS*}<5 






— 64 (all — fia) 



( 59 ) 
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with 


61 — - — -di, b 2 = 3 TuI/ 0 ’ b 3 = C 3 — 2^2) &4 = C* + 2 ^ 2 , 

65 = 6i5 pp ((u)) + b 2 R pp {p,p), b 6 = 2b 3 Spp{{u)) (60) 

and g having the same expression as in the 3-D algebraic equation. The two-dimensional polynomial 
solution of Eq. (59) is, as before, also written as 

a = £c A r\ (61) 

A 

Unlike the three-dimensional solution, however, the generators now consist of only five tensor 
groups: 

8 8 

T° = — — — , T 1 = S*, T 2 = S*ft-ftS*, T 3 = R*, T 4 = R*ft - hr* (62) 

o ^ 


for which there are five non-zero independent invariants, 

a 2 = {S* 2 }, w 2 = -{ft 2 }, {R* 2 }, {R*S*}, {S*B*ft}. (63) 

The fact that there are no other independent tensor generators [or invariants] can be verified using 
the 2x2 matrix identity: 


2 abc = bc{a} + a{bc} + ac{b} — b{ac} + ab{c} + c{ab} — 

c{a}{b} - a{b}{c} + ({ac}{b} - {acb})^ 2 *. (64) 


To obtain the solution to the algebraic equations for the anisotropy tensor, Eq. (59), a procedure 
similar to the one devised by Pope 20 is used. Three 5x5 matrices H*, 1 are defined: 


t. a _ i m 

l 3 2 


T^S* + S*T” - -{T^SM = 

3 A 

T"ft* - ft*T” = 


(65) 


for which A = 0 - 4. The elements of the matrices are determined from the above equations by 
making use of matrix relations stemming from the Cayley-Hamilton theorem. The 2-D tensor 
polynomial a = C X T X is introduced in both sides of Eq. (59). By making use of the above 
matrix identities the coefficients of the tensor polynomials are found to satisfy the following system 
of equations: 


C x = —gr 


Mia + Ms a + b 3 Y C n nl -b 4 Y C'JZ - Moa - b 6 £ C'F X 


( 66 ) 


The solution of this system of equations determines the model coefficients in the following algebraic 
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expression for the anisotropy tensor 


a = -2 C^t 

[C 2 S* + «2i + Qz)hgf 2 ra 2 S - S _ n S*)] 



—2 C's [R* + 6 ^/ir(R*fi - fi R*)] 

(67) 

suitable for two-dimensional mean flows. The eddy viscosities, C M and are given by: 



c _ hgfi/2 

^ 1 ~ §(&3<r) 2 cr 2 /i / 2 + 2(b 4 gf l T) 2 u 2 

( 68 ) 


r , _ b 2 gfi/2 

M 1 + 2{b 4 gfiT) 2 u 2 ' 

(69) 

The Q t coefficients 

are given by 



_ , b 2 r{s*R*> „ . . {s*R*mi 

Ql = 1 + u V + W r ‘- ? ] 

(70) 


0 - 1 1 2 ( ! >39 r ) V /l/2 . kk , 

Qi + 3\ + 2{b A ghT) 2 u^ Ql } 36j ^ 2 

(71) 


0 h l + 2{b A gfiT) 2 u 2 

3 M 3 2gf\To 2 

(72) 


with /i = (1 + b§gr / 6) 1 and f 2 — — begr/6) Note that the direct effect of compressibility as 

reflected in the baroclinic dyad occurs in Q\ and Q 2 - 

The high order of nonlinearity of the algebraic equations does not permit, in general, the 
construction of a fully explicit solution. Instead, an iterative approach is employed during the 
computations to generate the correct values. The algebraic solution is linearized by lagging the 
turbulence production term which contains the nonlinearity. 


C. Discussion 

To conclude this section some general statements regarding the behavior of the algebraic closure 
derived for the Reynolds stresses are highlighted. In Sec. Ill a closure for the Reynolds stress equa- 
tions was obtained. In the present section the fixed point solution of the modeled second moment 
equations, under the condition of structural equilibrium, was obtained. For two-dimensional mean 
flows the compressible algebraic stress model can be symbolically written as 

a = - I'to S* - u tl [S*Q -OS*] - ua[p - <S (2) ] - *& R* - v c tl [R*fi - Q R*]. (73) 

Note that products of the mean strain and the baroclinic dyad do not appear. This is due, for 
two-dimensional flows, to the relation T^S* + - |{T 17 S*}<S = -2{T y? S^}T°; the generator 

comprised of the mean strain and baroclinic dyad product is redundant. Examining the above 
expression, the following observations can be made: 

1. The first two terms are the same terms obtained in algebraic stress closure for two-dimensional 
[in the mean] incompressible flows. 
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2 . The first two eddy coefficients are functions of the relative strain and the relative rotation 
as is the case in incompressible flows. They are now also functions of the turbulent Mach 
number and the gradient Mach number. 

3 . Neither of these eddy coefficients depend on baroclinic effects. 

4 . It is seen that in the absence of mean velocity gradients that the turbulence is anisotropic 
due to the mean baroclinic dyad. This anisotropy manifests itself in the deviatoric as well as 
diagonal terms. 


It is useful to construct a simple example to see how the new effects influence the anisotropy. 
Consider a simple shear flow with a streamwise acceleration: let the Favre mean velocity, the mean 
density and the mean pressure gradients be represented by U\, i, Vp = />, 2, and VP = P, 1# 

2 11 

an = - g v i0 Uux + 2 U%1 ^1,2 + g Uiy2 p y 2 P,i (74) 

<*22 — 2 Vto U\,\ - - V t \ U12 + 2 Vt2 — V ft C/1,2 P,2 Pjl ( 75 ) 

1 2 

<*33 = j v t o C/1,1 — g ^<2 (76) 

<*i2 = - - ^vnU\,2 U\,\ - v*i 0 p,2 P,\ (77) 

Several observations regarding the above algebraic expression for the anisotropy can be made: 


1 . The expression is the first [that we know of] rigorous indication of the direct role the baroclinic 
dyad plays in determining the Reynolds stresses. 

2 . The expression indicates that, for arbitrary mean deformation, the mean baroclinic dyad 
contributes to the deviatoric portions of the Reynolds stress. 

3 . The expression also indicates that the baroclinic dyad also changes the relative magnitude of 
the normal stresses. This effect only occurs for mean deformations that are rotational. For 
an irrotational mean deformation the baroclinic dyad makes no contributions to the normal 
stresses. 

4 . For a uniform mean velocity the baroclinic dyad is a source contributing only to the deviatoric 
portions of the anisotropy tensor. 

5 . The expression indicates the inapplicability of any heuristic gradient transfer arguments for 
the Reynolds stresses in flows with important gradients of mean density and pressure. 


While these results indicate the inapplicability of any form of eddy viscosity model for the stresses 
in compressible flows with arbitrary large density and pressure gradients some qualifications are in 
order. The presence of the baroclinic dyad is likely to be important only in rapidly accelerating 
aerodynamic flows or in combusting flows where one can expect the mass fluxes to be important. 
In the absence of these effects it appears that the parameterization of the Reynolds stresses in 
terms of powers of the mean deformation with modifications according to the compressibility of the 
fluctuations as those indicated in Sec. Ill is appropriate. 
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V, Computational investigations of free shear flows 


The theory and results presented in the previous sections are now implemented over a very wide 
range of mean flow Mach numbers in the simulation of free shear flows. Simulations using second- 
order moment (SOM) closures as well as the quasi-explicit algebraic models are conducted. The 
numerical experiments are constructed with the intention of investigating several different issues of 
relevance to the prediction of compressible turbulent flows for engineering purposes. 

In addition to assessing the sensitivity of the pressure-strain model to the different pressure- 
dilatation models we compare the computational results to what is expected from laboratory and 
numerical experiments. Of particular interest is the well known effect of compressibility on reduction 
of the spread rate of the mixing layer. The main objective is to assess the effectiveness of the 
compressibility corrections in reproducing the reduced growth rate of the free-shear layers with 
increasing Mach number - a phenomenon which is well documented experimentally. 44 We also 
study the adequacy of the algebraic stress closure with the results of the SOM simulations. In 
particular the Reynolds stresses and the anisotropy computed using both of these procedures are 
compared to solutions without the compressible corrections in order to assess the effects of the new 
pressure-strain closure on the anisotropy. Our intention is to further investigate the notion that 
the dramatic reduction of the mixing layer growth rate is due to the effects of compressibility on 
the pressure-strain and, consequently, its effect on the reduction of the turbulence shear stress. 


A. Numerical method 

The simulations are conducted using a finite difference scheme 45 second-order accurate in time 
and fourth-order accurate in space. The governing equations are integrated explicitly in time using 
the predictor-corrector finite difference scheme developed by Gottlieb and Turkel. 46 The Gottlieb- 
Turkel scheme is a higher order accurate variant of the MacCormack 47 predictor-corrector method. 
During a numerical sweep, the inviscid fluxes are alternatively differenced backward in the predictor 
step and forward in the corrector step. Fourth-order central differences are used for the viscous 
and heat flux terms as well as for the derivatives in the source vector. To maintain stability, it 
is required that the CFL number be less than 2/3. To prevent numerical oscillations in regions 
of large gradients a smoothing scheme devised by MacCormack and Baldwin 48 is employed. The 
method, as outlined in Ref. 45, adds artificial viscosity that is very small everywhere except in the 
regions where the pressure oscillates. 

All simulations are conducted on a uniformly spaced grid in the computational domain. By 
means of a coordinate transformation the mesh is transversally compressed in the physical domain 
in the region corresponding to the mixing layer. 

The physical domain is a rectangular box defined by the set of points (x,y), in which x rep- 
resents the streamwise coordinate and y the transversal coordinate. The grid overlaying the com- 
putational domain of size x8 w x y8 w has 100 x 60 points, where the vorticity thickness & w = 
(u\ — u 2 ) / (d(u) / dy) max . In this nondimensionalization, the reference length scale is the magnitude 
of the vorticity thickness at the inlet. The initial profile for the mean axial velocity is adjusted such 
that the resulting nondimensional vorticity thickness at the inlet is equal to one. The values of the 
physical dimensions are (x,y) = [80,20] for the mixing layer. To evaluate the grid sensitivity the 
number of grid points is increased by a factor of 1.5. The change in the steady-state values of the 
peak of the turbulent stresses is less than 2 percent. 
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To accelerate convergence to the stationary state, a local time stepping technique is used. The 
convergence criterion is imposed so that the global averaged residual profile attained is station- 
ary for each dependent variable. Although more stringent criteria can be sought, it is known, 49 
that predictor-corrector schemes are limited in their ability to achieve very high rates of residual 
reduction. 

Due to the nonlinearity of the destruction terms, the k — e equations display a stiffness which can 
either generate numerical instabilities or increase the computational time. In order to avoid these 
inconveniences the turbulence source terms are treated implicitly. The Ar — e destruction terms are 
decoupled by suitable manipulation of the e/k ratio, treated as a known quantity from the previous 
time step. 


Initial and boundary conditions The initial fields are obtained by propagating the inflow 
conditions throughout the entire domain; hence, the flow has to sweep the domain at least one 
time to obtain meaningful results. For ease of computation, the inflow initial conditions (IC) for 
the flow variables are assumed to be a smoothed step (hyperbolic tangents for the axial velocity or 
species) or a smoothed hat profile (for turbulence quantities). Uniform profiles are assigned at the 
inlet to the static pressure, to the static temperature and, by virtue of the equation of state, to the 
static density. 

The boundary conditions (BC) are set according to the elliptic nature of the problem on all 
four boundaries. The inflow BC are fixed (Dirichlet) for all primary variables in the supersonic 
and subsonic regions with one exception. For the portions where the flow is subsonic the pressure 
is allowed to adjust to the characteristic waves through a Neumann boundary condition. At the 
outflow and outer boundaries zero gradient (Neumann) conditions are applied due to their non- 
reflective properties in relation with the outgoing waves. In the mixing layer, the static pressure, 
temperature and density in the two free stream layers are matched. This isolates the contributions 
to the layer growth to those due solely to the variation of the reference Mach number. 


Free-shear flow parameters The magnitude of the effects of compressibility are often param- 
eterized by the convective Mach number. 50 In our formulation the expression for the convective 
Mach number is M c = Af( 1 - r v )/ 2, where r v = U2/U1 is the axial velocity ( u ) ratio. A wide range 
of values of M c including both subsonic and supersonic regimes are considered: 0.2 < M c < 2.0. 
The convective Mach number is varied by keeping the velocity ratio constant, r v = 1/2, and varying 
the reference Mach number M . The reference Reynolds number is Re = 5 X I0 6 , while all other 
nondimensional parameters are kept constant. The subscript 1 refers to the high-speed stream 
value, while the subscript 2 refers to the low-speed stream value in the mixing layer. In the results 
presented, the spatial coordinate, for the hydrodynamic variables, is given by 77 = £ y ^ x °' 5 \ where 
u = . Only the mixing layer simulations are presented in this article. Simulations for the jet 

configuration can be found in Ref. 37. 

After an initial period of flow development a linear growth rate of the shear layers is attained. 
In the fully developed regime, when a linear growth rate has been established, the mean velocity 
and the normalized Reynolds stresses display self-similar behavior. The spread rate of a turbulent 
shear layer, in the self-similar region, is conventionally expressed as dSJdx = C s {\ - r„)/(l + r v ) 
where S u could be either the 10 percent visual thickness of the shear layer based on the normalized 
velocity profile or the vorticity thickness, and Cs is a constant (approximately). As the present 
article is concerned with compressible corrections to the turbulence moment equations no plots of 
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the mean flow quantities are given; the mean field variables are quite typical and can be seen in 
the literature. 11,51,52 


Model calculations Computations are performed with several different forms of closures. The 
presentation of the results is much clearer if we designate precisely what equations are used for 
what calculations. The central methodology is based on the models of Ristorcelli . 8 Thus KE-Rist 
represents the k - e model equations Eqs. (4, 5, 6, 7, 18, 20), with the Ristorcelli compressible 
closures Eqs. (34, 45). KE-Sarkar denotes the k — e model equations, with the Sarkar compressible 
closures Eqs. (33, 43). KE-Zeman is the k-e model equations, with the Zeman compressible closure 
Eq. (44). ARSM denotes the k — e model equations Eqs. (4, 5, 6, 7, 18, 20), with the 2-D ARSM 
compressible closures Eqs. (38, 48, 60, 67-72). ARSM-Sarkar denotes the k — e model equations, 
with the 2-D ARSM compressible closures Eqs. (36, 48, 60, 67-72). SOM represents the second 
order model equations Eqs. (4, 5, 6, 7, 15, 20), with the compressible closures Eqs. (28, 38, 42). 

B. Mixing layer simulations 

The primary concerns of this article are with 1) a representation of the effects of compressibility 
as they appear in the pressure-strain covariance and 2) the construction of an algebraic Reynolds 
stress model useful for engineering calculations. Studies related to both these issues, for the M c = 
1.07 mixing layer configuration, are shown in the first two figures. The Reynolds stresses and the 
production/dissipation ratio obtained for ARSM and SOM calculations are shown in Figs. 1 and 
2. Also shown, as a double check, are curves based on an a priori evaluation of the algebraic 
stress model; the data from a SOM calculation is used as input to the ARSM. The neglect of 
turbulent transport in the algebraic stress model is reflected in differences in the normal stresses at 
the centerline. In general the a priori evaluated ARSM curves follow closely the computed ARSM 
values indicating that for this flow the neglect of transport in the moment transport equations is 
justified. The algebraic approximation is, as has been verified at several different convective Mach 
numbers, suitable for an investigation of trends in this flow configuration. 

Also shown in Figs. 1 and 2 is the SOM computation without compressibility corrections [NC - no 
compressibility], i.e. d\ and ^2 are set to zero. The effect of the compressibility corrections manifest 
themselves as decreases in the centerline values of {uu),{uv) and (uu) of approximately 16%, 25% 
and 24% respectively. This is consistent with that seen in DNS. In addition, the streamwise normal 
stress suffers a smaller reduction than the other Reynolds stresses. This, as will be seen, will 
manifest itself in an increase in the streamwise normal anisotropy. 

The production/dissipation ratio is shown in the bottom of Fig. 2. The average level of P/e for 
the compressible case is smaller than the incompressible case; this is consistent with the reduced 
turbulent kinetic energy of the compressible flow. While the overall level of the production is 
smaller it is interesting to note that at one location in the mixing layer, the centerline, P/e is 
nominally higher. It has been conjectured that this might be due to the higher relative strain rates 
compressible flows may be able to sustain. This is an issue that further DNS may resolve. 

Figure 3 displays the vorticity thickness. Vorticity thickness spread rates as predicted by the 
SOM simulations, with and without (NC) the compressibility correction, are shown on the top of 
Fig. 3. It is seen, that even with such modest contributions from compressibility, scaling as they do 
with M 2 , a sizable change in growth rate is seen. The variation of M t with M c is shown in Fig. 5. 
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Figure 3 also displays the different components of the pressure-strain tensor computed from the 
SOM simulation with and without (NC) the compressibility correction. The different components 
of the pressure-strain are all reduced by about the same amount. This is precisely the behavior 
seen, in both direction and magnitude, in the recent DNS 15 - see figure 4.18 from Ref. 15. 


Reynolds stress anisotropy In Fig. 4 the computed values of the centerline Reynolds stress 
anisotropies are shown. At low values of compressibility the anisotropy has values similar to that of 
simple incompressible shears. At higher levels the anisotropy behaves very much like the compress- 
ible shear DNS; more of the energy of the turbulence is in the streamwise component as made clear 
by the increase in an. This is at the expense of the crosstream component, (tw), which is reduced. 
This is, as shall be discussed further, due to the reduction of intercomponent energy transfer due 
to the compressible aspect of the pressure-strain. There is also a reduction in the shear stress as is 
seen in compressible DNS. 


Comparison with laboratory Reynolds stresses The effects of the new compressible models 
on the anisotropy is shown in Fig. 4. Unfortunately, in the laboratory situation, the spanwise 
Reynolds stresses are rarely measured. As a consequence the anisotropy is not known and one 
returns to primitive variables. 

The maximum values of diverse Reynolds stresses in the mixing layer are plotted as function of 
the convective Mach number in Fig. 4. In this figure a u and a v represent, respectively, the (au) and 
(uu) Reynolds stresses nondimensionalized by the square of the velocity difference across the layer. 
Also shown are the experimental results of Elliot and Samimy. 53 Underlying this comparison is, of 
course, the implicit assumption of self-similarity. The maximum Reynolds stresses decrease with 
compressibility. In the present simulations, however, the peak values of the turbulence intensities 
(axial and transverse) do not vary as strongly as in the laboratory findings of Goebel and Dutton. 54 
In fact the present simulations are closer to the results of either Ref. 54 or Ref. 53 than either Ref. 54 
or Ref. 53 are to each other . 

There appears to be a discrepancy. A survey of the diverse compressible DNS indicates that 
an increases with compressibility. This implies that (uu) decreases more slowly than (uv). [Figure 
9 of Ref. 14 is the most comprehensive presentation of this trend.] This does not appear to be the 
case in all laboratory experiments. In the mixing layer of Ref. 53 (uu) decreases more rapidly than 
(uu) for all M c implying that a xx decreases with increasing compressibility contrary to DNS results. 
In Ref. 54 (uu) decreases more rapidly than (uu) for low Af c ; while at moderate Af c , (uu) decrease 
less rapidly than (uu) in accord with DNS results. No explanation for this discrepancy between 
laboratory and numerical data is known. 


The pressure-strain tensor The effects of compressibility on the Reynolds stresses are believed 
to be the source of the dramatic influence on the spread rate as the convective Mach number 
increases. The physical mechanisms responsible for these phenomena were once believed to be the 
compressible dissipation and pressure-dilatation which caused a decrease in the Reynolds stresses 
due to an overall decrease in k . This was accompanied by no substantial changes in the anisotropy. 
As discussed in Sec. I, recent simulations suggest that compressibility causes a turbulent shear stress 
reduction due to a reduction in its associated anisotropy. This is a statement that compressibility 
manifests itself structurally, not energetically. The reduction in the anisotropy leads, of course, to 
a reduction in k and ultimately in the mixing layer growth rate. 
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More recent simulations, 9,14-16 all seem to indicate that the major changes in the shear stress 
are due to the effect of compressibility on the pressu re-strain covariance. Figure 6 is an illustration 
of this idea: Fig. 6 depicts the different components of pressure-strain tensor and its trace [the 
pressure dilatation] at the centerline of the mixing layer. Figure 6 is patterned after Vreman et 
ai ’s 14 figures 3 and 9. The behavior of all the components of the pressure-strain, shown in Fig. 6, 
are very much in accord with the figures of Ref. 14. The only notable difference is the relative levels 
of 1122 and II33. The trends are very much the same and the relative behavior - a commensurate 
across the board reduction with M c - is very much in accord with Ref. 15. 

The behavior of the pressure-strain illustrated in Fig. 6 is very much in keeping with the current 
understanding of the mechanism for the reduction of the turbulence energy by the reduction of the 
shear stress: the only source for the kinetic energy is in the (uu) component and its production 
is proportional to ( uv ). The reduction [in magnitude] of II n and II 22 means that less energy 
generated in the (uu) equation is transferred to (vv) [by 1122 ]* As (uv) has no production (uu) is 
much smaller. As consequence of the fact that the production of (uv) is proportional (vv) much 
less (uv) is produced and much less (uu) is produced. This can be seen by the considering the 
truncated forms of the second-order moment equations: 

~ - (uv)Ui ,2 + Iln + .... (78) 

§- t ( vv ) ~ n 22 +.... (79) 

-^(uu) ~ - (vv)Ui ,2 + II 12 + (80) 

The thesis of Freund 15 provides a schematic of the above process. The behavior of the Reynolds 
stresses and the pressure-strain indicated by Figs. 2 and 6 is consistent with the mechanism for the 
reduction just delineated. Specifically, the fact that less of the (uu) energy is transferred to the 
other components of the Reynolds stress means that the anisotropy, an = 2 &n, must increase as 
seen in all the DNS. Note also, as might be expected for a quasi-equilibrium flow, that the pressure- 
dilatation as derived by the asymptotic procedure in Ref. 8 is very small; negligible in comparison 
to the pressure-strain. Yet its use to obtain the deviatoric portions of the pressure-strain produces 
a very sizable reduction in the shear stress and the growth rate. 


Gradient Mach number Sarkar 16 and earlier work cited therein have drawn attention to the 
gradient Mach number as a potentially useful parameterization of the effects of compressibility. 
Many of the effects of compressibility, as indicated by diverse DNS, have been observed to become 
more apparent as the gradient Mach number increases. Shown in the bottom of Fig. 5 is the 
gradient Mach number using the definition M g — \M t a(k)/e a . The definition of the gradient Mach 
reflects the fact that the Kolmogorov scaling, e ~ (k) 3 ! 2 ft has been used to eliminate the length 
scale in the definition M g — oijc. In DNS, the length scale t is determined from the two-point 
correlation; no such opportunity occurs in single-point closures. Both the DNS, for example figure 
4.29 in Ref. 15, and the SOM simulation indicate a similar decrease in rate of increase of M g with 
M c . The decrease seems to occur at lower M c in the DNS. This is likely to be related to the different 
length scales used. 

One of the noteworthy observations made of the effects of compressibility is their tendency to 
saturate. 11 In the bottom of Fig. 5 the rate of increase of the maximum turbulent Mach decreases 
at higher convective Mach numbers. This is in line with other observations. 11 


26 



Kolmogorov scaling While in DNS a length can be determined from the two-point correlation, 
no such possibility exists for single- point closures. Yet a length scale, reflecting the two-point nature 
of the turbulence problem, as seen analytically in Ref. 8, is required. The length scale appearing 
in the Kolmogorov inertial range scaling, e ~ (k) 3 ^ 2 /£ } was used in Ref. 8. 

There is a proportionality coefficient, a, in the Kolmogorov scaling a = I s £/(2k/3) 3 ^ 2 . For infi- 
nite Reynolds number isotropic turbulence a ~ 1. For finite but large Reynolds number anisotropic 
turbulence undergoing deformation the Kolmogorov scaling is likely to be useful but a is likely to be 
a flow dependent quantity. Sreenivasan 55 has made some studies of the variation of the Kolmogorov 
scaling coefficient: he has found values of a ranging between 0.4 — 2 for diverse incompressible sim- 
ple shear flows. To allow for the expected variability of a two different values of are used in many 
of the simulations. The value of a does not change any of the trends and does not seem to have 
too strong an effect on the energy at the centerline (see Figs. 4 and 5). It does however effect the 
mixing layer spread rate as can be seen in Fig. 8. 


Mixing layer growth rates Settles and Dodson 44 have compiled a very large number of exper- 
imental results. 50 ’ 54 ’ 56 " 61 These are shown in Figs. 7 and 8. The scatter in the data reflects the fact 
that different experiments are done in different wind tunnels, with different inlet configurations and 
with different reservoirs. The computational curves shown in Figs. 7 and 8 reflect two different but 
complementary investigations. As the “Langley curve” 62 has been a popular benchmark it has also 
been included. The well known reduction in the mixing rate with increasing M c is seen in all the 
data. 

Figure 7 reflects computations using the “energetic” approaches. A k — e scheme has been used 
to calculate the mixing layer growth and the only compressible corrections are in the k equation; 
the Reynolds stresses are closed with the usual incompressible eddy viscosity [Boussinesq] approxi- 
mation. Figure 7 indicates that the results based on the KE-Zeman and KE-Sarkar schemes capture 
the mixing layer reduction. The reduction of the mixing layer growth rate has been captured by 
similar mechanisms in both models; most of the suppression is provided by a substantial amount of 
additional dissipation that comes from the compressible dissipation models. As was pointed out in 
Sec. I, the compressible dissipation, as indicated by several DNS and by asymptotic analysis, 8 does 
not play an important role in these classes of flows. Any arbitrarily dissipative term added to the 
k equation, appropriately calibrated and scaled, can produce very good agreement for the mixing 
layer growth. Moreover, as with all energetic approaches, there is no possibility of accounting for 
the very important structural changes that appear in the anisotropy. 

As indicated by Fig. 7 the pressure-dilatation model of Ristorcelli 8 by itself, with a = 2, accounts 
for a nominal suppression of the growth rate. The compressible dissipation being negligible in 
analysis of Ristorcelli 8 is not included in this calculation. From incompressible DNS, 63 it can be 
argued that a value of a % 1 might be more appropriate. As the pressure dilatation scales with a 2 
a value of a % 1 would decrease the pressure-dilatation effects by a factor of four. 

Figure 8 reflects a computation using the compressible algebraic Reynolds stress model. This is 
exactly the same computation as given in the Fig. 7 with the exception that the algebraic Reynolds 
stress approximation now includes the effects of compressibility in the pressure-strain covariance. 
Note that there are substantial changes on the layer growth rate prediction. The change is more 
drastic for the ARSM modeling than the ARSM-Sarkar modeling; this is to be expected as the 
majority of the growth rate reduction using the earlier Sarkar modeling was built into a dissipative 
term (which does not contribute to the modeling of the pressure strain) and the relative change is 
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small. 


The algebraic model built upon dilatational closures of Ristorcelli 8 shows improvement over 
the results from a simple k — e scheme. Figure 8 depicts the outcomes of the computations with 
the Ristorcelli 8 based algebraic closures for two values of the Kolmogorov scaling coefficient. The 
results for a = 0.4 with both the pressure-strain and pressure-dilatation provide a modest decrease 
in the spread rate - matching that predicted with simply the pressure-dilatation with a = 2. The 
ARSM calculation predictions are substantially reduced when a = 2. If a were an undefined ad 
hoc constant one might be tempted to set it so as to match the reduction in spread rate indicated 
by the data. Given that little is known about a [though it is sure to be in the range 0.4-2. 0 and 
most likely a « 1] these issues must be explained by DNS. 

It is concluded that the present pressure-strain modeling method based on the extension of 
the well-established incompressible procedure, in which the trace-free constraint is relaxed, does 
rationally account for a significant portion of the reduction in shear stress and growth rate. This 
procedure is considered a leading order contribution to the compressible turbulence shear stress 
problem. The possibility that additional compressible corrections to the pressure-strain, addressing 
physics not able to be accounted for using this constitutive relation methodology, may require 
consideration. 


VI. Some issues in compressible turbulence modeling 

In our effort to obtain a closure for the effects of compressibility a few issues, already alluded 
to, have become clearer. These issues are now highlighted. 


The baroclinic dyad The procedure invoked to obtain the algebraic Reynolds model indicates 
that the baroclinic dyad is a source of turbulence. This, of course, can be anticipated by inspection 
of the second-order moment equations. The consequence is that an eddy viscosity representation 
for the Reynolds stresses is, from first principles, inappropriate for classes of flows in which the 
mass fluxes are important. 

The presence of the baroclinic dyad is likely to be important only in rapidly accelerating aero- 
dynamic flows such as through shocks or in hypersonic situations. The baroclinic dyad is also 
likely to be important in combusting flows where one can expect the mass fluxes to be important. 
In noncombusting supersonic flows it appears that a parameterization of the Reynolds stresses in 
terms of the mean deformation with modifications for the compressibility of the fluctuations as 
derived in Sec. III. A appears appropriate. 


A length scale for single-point models of compressible turbulence In Ref. 8 the effects of 
the non-zero divergence of compressible turbulence was parameterized by several two-point integrals 
made nondimensional by a length scale. For single-point turbulence closures the length scale chosen 
is typically that appearing in the Kolmogorov scaling, € ~ k 3 / 2 /£, which is often taken to be the 
integral length scale of the longitudinal correlation. 

For compressible turbulence of interest to supersonic aerodynamic flows the cascade mechanism 
will be comprised of the usual nonlinear solenoidal modal interactions. The Kolmogorov scaling 
is expected to be valid in the weakly compressible limit. However, given the observed effect of 
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compressibility on the length scale seen in simulations (see figure 4.28 of the recent report 15 ) one 
might expect the coefficient of proportionality, <*, to be a function of compressibility. That this 
might be the case is also suggested by the appearance of different instability modes in compressible 
flows. Studies similar to the incompressible studies shown in Ref. 63 appear to be both interesting 
and very relevant to issues of compressible turbulence modeling. 

It should be understood that the Kolmogorov scaling coefficient, a = e s £/(2k/3) 3 / 2 , is not a free 
parameter in the sense used in many turbulence model closures. The Kolmogorov constant appears 
in definite physical relationship having a precise mathematical definition, a clear interpretation in 
terms of flow physics, and substantial theoretical and numerical substantiation. It is, in principle, 
measurable for any flow. This in contrast to an ad hoc modeling constant; such a constant is often 
determined by calibrating the results of the model equations to some experimental data. Such a 
“calibration constant” has neither a precise mathematical definition or a clear interpretation in 
terms of flow physics and is dependent on the set of model equations used to compute the flow. 


Experimental differences in numerical and laboratory data There appears to be a unani- 
mous agreement in the DNS that, with increasing gradient or convective Mach number, an increases 
while 022 and ai 2 both decrease. This trend does not appear in the laboratory mixing layer exper- 
iments. At low M C} (uu) appears to decrease more rapidly, with increasing compressibility, than 
(vv). This implies that an decreases with compressibility in contradistinction to DNS results. It is 
crucial to understand the cause for the differences between the DNS and the laboratory flows of the 
{uu) Reynolds stress behavior; turbulence models are based on intuition gleaned from DNS data 
but are used to predict engineering flows that are more similar to laboratory flows. The present 
computational results are closer to either of the laboratory experiments of Ref. 53 or Ref. 54 than 
the laboratory experiments of Ref, 53 or Ref. 54 are to each other. Earnest speculation on the 
source of the differences in the two experiments and facilities is required. 

VII. Summary and conclusions 

Progress towards the development of a compressible turbulence closure, starting at the level 
of second-order moment equations, has been described. Modeling from the second-order level 
accommodates important structural changes that appear in the anisotropy and are a feature and a 
function of compressibility. In the second-order moment equations the compressible contributions to 
the pressure-strain covariance have been obtained. The pressure-strain has been closed by assuming 
that, as is consistent with the weakly compressible limit, it can be modeled as a tensor polynomial 
linear in the Reynolds stresses. The difference from the incompressible case is that the trace of 
the compressible strain is not zero; it is set equal to the pressure-dilatation for which models exist. 
The compressible pressure-strain closure features a dependence on several turbulence descriptors: 
the turbulent Mach number, the relative strain, the gradient Mach number, the production and 
the dissipation. As a consequence the coefficients in the compressible pressure-strain closure are 
not constants but functions of parameters the turbulence and its compressibility. 

In addition to devising a closure for the pressure-strain, a closure for the mean acceleration/mass 
flux terms appearing in the Reynolds stress equations has also been developed. For the flows 
studied in this article, limited as they are to the flows for which experimental data are available, 
the acceleration/mass flux moments are not important. The mass flux terms will be important in 
the combusting or hypersonic flows which motivated the thesis. 37 Further development of the mean 
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acceleration terms requires experimental data for this class of flows. 

Having closed the compressible Reynolds stress equations standard tensor representation theory 
has been used to produce a compressible algebraic Reynolds stress model useful for flows near 
structural equilibrium. A noteworthy feature of this portion of the work is that the role that the 
mean pressure and mean density gradients plays on the Reynolds stresses is immediately seen. As 
a consequence, it is seen that the baroclinic dyad can, in the absence of mean velocity gradients, 
contribute to the anisotropy of the turbulence. It is also seen that the mean bulk dilatation 
contributes to the anisotropy in a way that is quite different from an irrotational mean strain. 

The mathematical results developed here have been implemented in mixing layer computations 
spanning a wide range of mean flow Mach numbers. In this article the discussion has been limited 
to the compressible mixing layer for which a sizable amount of literature exists. The calculations 
presented have been organized along two themes: 1) an investigation of the effects of compressibility 
as related to the compressible pressure-strain and the Reynolds stresses and 2) a validation of the 
algebraic Reynolds stress model predictions. 

The computations with the compressible pressure-strain indicate that the present modeling 
[which does not have undefined tunable constants] produces precisely the behavior seen in the 
DNS 14,15 ; there is a commensurate reduction, with increasing convective Mach number, of all the 
components of the compressible pressure-strain tensor. The changes in the pressure-strain lead, as 
is established by DNS of several different flows, 14 * 16 to changes in the anisotropy of the Reynolds 
stresses. The changes predicted by the modeling for the normal and shear anisotropies are very 
consistent, in trend, with DNS data and especially with the DNS of the mixing layer 14 — the flow 
configuration most similar to the one treated in this article. Comparison between the computa- 
tional predictions and laboratory results for the Reynolds stress and their behavior with increasing 
compressibility is not as good as might be hoped. This is due to unknown experimental issues: the 
agreement between even the different laboratory experiments is poorer than the agreement between 
the numerical and experimental results. 

All the laboratory data do however agree on the trend, with increasing compressibility, of the 
mixing layer growth rate and its kinetic energy: it decreases. The computational experiments 
conducted indicate that sizable reductions in the mixing layer growth rate accompany changes in 
the anisotropy of the turbulence due to the compressible aspects of the pressure-strain covariance. 
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Appendix A 


The Poisson equation for the fluctuating pressure is obtained by taking the divergence of the 
Navier-Stokes equations. Both instantaneous and averaged forms of these equations are necessary 
in the derivation. The subtraction of the mean equations from the instantaneous equations and 
simple term manipulations involving the continuity relation provide the desired Poisson equation 5 : 
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The solution for the fluctuating pressure can be determined by using Green’s functions or Fourier 
transforms. With respect to the incompressible counterpart, the compressible Poisson equation is 
severely complicated. An equivalent form which resembles more to the incompressible equation can 
be obtained as: 
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The last term on the right hand side (RHS) converts, when the fluctuating pressure solution is used 
to determine the pressure-strain correlation, into a viscous interaction tensor whose trace is equal 
to the dilatational dissipation. In the absence of walls this term can be omitted in high Reynolds 
number flows. The acoustic term (the first term on the RHS) is difficult to be taken into account 
in the present analysis. If it is assumed that the pressure fluctuations are caused by turbulence 
only then this term is negligible. The condition p'/p 1 allows the neglect of the second term 
on the RHS. The remaining terms are the return-to-isotropy and the rapid part for which known 
modeling principles can be applied in the limit of low convective Mach numbers. 


Appendix B 

It well known that for linear pressure strain forms it is impossible to satisfy realizability con- 
ditions - the requirement that the eigenvalues of the Reynolds stress tensor remain positive. Sat- 
isfying realizability is a very practical computationally stabilizing requirement. Our experience 
with computations in complex flows indicates that realizability is very useful. The model is now 
made (weakly) realizable following methods suggested by Schumann 64 and Lumley 34 and detailed 
by Shih and Shabbir. 65 Let F — 1 + 27/77/8 + 9/7/4 is a parameter involving the second invariant 
77 = dijdji and third invariant III = of the Reynolds stress anisotropy tensor. 

Then the following asymptotic behavior for the pressure strain-model ensures that realizability is 
satisfied: 



= CF a 


as F — > 0 
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(83) 


diup}j 

dx q Xpe9e 


-> 0 


as F — *■ 0 


where the index e indicates that the relations are written in the principal axes of (u"u'j). The 
computational form of the model with the additional parameters necessary to enforce realizability 
so necessary for computational stability is 


n* — p (*j — —C\p ta i: A r F ar 


m [(| + f *) 

[1 - C 3 + 2 ,d 2 ] [a,„5* «u» + S* p ((u))a p j - |s^«ti»a F ,«y] - 
[1 ~ C A - 2 d 2 ][a ip n pj ((u)) - Q ip ((u))a pj ] + ^d 2 5 pp ((r t »a, i ] B r F 0r 


(84) 


with Cz = (5 — 9C2)/11 and C 4 = (l-f7C2)/ll. The value for the constant C? will be the same as in 
the incompressible model to preserve consistency in the zero Mach number limit, that is C 2 = 0.45, 
The parameters are a r = 0.1, 0 r = 0.5, A r = min (F _Qfr , 0.1~ Qr ) and B r = min(F“ /3r , 0.1“^ r ). 

Using this modified form of the pressure-strain model the ARSM coefficients become: 61 = 
i - B r F^(i + §d,), 6 2 = 63 = 1 - B r F^(l -C 3 + 2 d 2 ), 64 = 1 - B r F^(i - C 4 - 2 d 2 ), 

and 


9 = 


Ar F a 'C^ + C £2 - 2 + (2 - + ~^+ 

f a p tt a Dt 


2 T 

3 


(l-2d 2 B r F^)S pp ((u)) + 2 


M + j/d - p f e 

pi* 


-l 


(85) 
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Figures 


FIG. 1. Comparisons between ARSM and SOM calculations for the mixing layer, M c = 1.07; (a) 
streamwise normal Reynolds stress; (b) Reynolds shear stress, [a = 1.2] 

FIG. 2. Comparisons between ARSM and SOM calculations for the mixing layer, M c = 1.07; (a) 
crosstream normal Reynolds stress; (b) turbulent kinetic energy production over dissipation. 
[a= 1.2] 

FIG. 3. Influence of the compressibility correction in SOM calculations for the mixing layer, M c = 
1.07; (a) the vorticity thickness; (b) components of the pressure-strain tensor, [a = 1.2] 

FIG. 4. Variation with convective Mach number in the mixing layer; (a) centerline Reynolds stress 
anisotropies; (b) centerline Reynolds stresses, [a = 1.2] 

FIG. 5. Variation with convective Mach number in the mixing layer; (a) maximum turbulent Mach 
number; (b) maximum gradient Mach number. 

FIG. 6. Variation with convective Mach number for the mixing layer; the pressure-strain tensor. 

[a=1.2] 

FIG. 7. Effect of compressibility on normalized thickness growth rate, k — e calculations of a mixing 
layer, [a = 2] 

FIG. 8. Effect of compressibility on normalized thickness growth rate, ARSM calculations of a 
mixing layer. 
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